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Appendix I 

TGP - TRAJECTORY GENERATOR PROGRAM 

I. 1 INTRODUCTION 

The Trajectory Generator Program (TGP) is a digital computer 
program written in FORTRAN IV to run on IBM 7090/7094 computers. The 
program is used in mission analysis and planning to generate a preliminary 
trajectory {or trajectories) from launch to orbit injection followed by several 
options for performing orbit transfer with or without rendezvous. 

The program employs analytical equations for determining launch 
conditions and the atmospheric phase of ascent. The vacuum phase of ascent 
(to orbit injection or intercept) is generated with an explicit guidance algorithm. 
All orbital burns are considered as impulses and are planned by two different 
techniques to provide alternate solutions. The program is based on a spherical 
earth model which allows for the rapid generation of coast segments with Kepler 
type equations. 

This program was written to demonstrate feasibility of the Prelim¬ 
inary Trajectory Generation Mode of the On-Board Mission Planning Function 
(see Section V of the main report). TGP is used to generate initial trajectories 
for the Trajectory Optimizer Program (TOP) which is described in Appendix 

II. 

1.2 PURPOSE AND USE 

TGP was written to demonstrate preliminary feasibility of on-board 
trajectory generation for the QRGT concept. The program also demonstrates 
the capability of the generalized guidance algorithms, developed during this 
study, for mission planning purposes. TGP represents the "Explicit" aspect 
of the "Optimal-Explicit" approach to on-board trajectory planning which is 
outlined in Par. 5.4 of the main report. 

. ■; This program is employed to generate preliminary trajectories for 
r various missions* that may employ the QRGT concept. These preliminary tra¬ 
jectories are specified a set of parameters which identify: launch conditions, 
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ascent to orbit, and all major orbit burns. These parameter sets then become 
the input to the Trajectory Optimizer Program-TOP (Appendix II) where they 
are systematically varied in order to minimize some mission performance 
function (e.g. total AV or total mission time). TOP then represents the "Op¬ 
timal" aspect of the "Optimal-Explicit" approach to trajectory planning. 

1.3 PROGRAM FEATURES 

The present version of TGP makes use of many subroutines gen¬ 
erated in different subtasks of Task II during this phase of QRGT. This has 
resulted in some duplication and program inefficiency. These deficiencies can 
be corrected when a more general and sophisticated program (which will be 
necessary for detailed feasibility studies) is written during Phase II. 

The present program employs fairly simple planning logic and 
launch routines in conjunction with quite sophisticated ascent and orbit maneu¬ 
ver routines to provide the following features: 

• Direct ascent, ascent to orbit and orbit maneuvers (with or 
without rendezvous) can be generated. 

# • Spherical gravity model. 

• Kepler arcs for coast phases. 

• Launch conditions options for minimum out-of-plane angle 
or for a specified launch time. 

• Maximum and minimum launch azimuth constraints included. 

• Ascent to orbit or intercept with an optimal explicit guidance 
algorithm (UP 1 or UP 2 Subroutines, which are really two 
entries to a slightly modified version of OP-EX, the algorithm 
used for powered-flight guidance and described in Appendix 
III). . 

• Two-burn (RENDI Subroutine) or three-burn (YPG Subrou¬ 
tine) orbit transfers possible. 

• Orbital burns treated as impulses. 

• Mission can be planned from launch or from a parking orbit. 

• The two burn rendezvous routine (RENDI) employs a search 
routine to locate local minimum- AV transfers. 
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1.4 PROGRAM OPERATION 

The over-all operation of TGP is illustrated in Figure 1-1. The 
program consists of an executive routine MAIN and seven major subroutines: 
LAUNCH, UP 1, UP 2, EPHEM, RVVR, RENDI and VPG. The MAIN program 
reads in all data, performs necessary conversions, sets up program paths, 
initializes and calls subroutines, and writes all output data. The LAUNCH 
subroutine is used to determine the launch time and launch conditions. The 
ascent to orbit or intercept phases is calculated in UP 1 and UP 2 respec¬ 
tively* Ephemeris for the spacecraft and target is updated with EPHEM. The 
RVVR routine is used to compute the required velocity necessary for an orbital 
burn. Two-burn orbital rendezvous maneuvers are generated with RENDI 
while VPG generates several different rendezvous trajectories with three (or 
more) burns per trajectory. 

For each successful run, the printout includes all parameters nec¬ 
essary to specify each trajectory. These include: launch time and azimuth; 
ascent cut-off state, time, and AV remaining; position, velocity, time and 
A V at each orbital burn; and other related data. 

1.4.1 MAIN PROGRAM 

The executive routine MAIN represents the planning logic and pre¬ 
scriptions of the Preliminary Trajectory Generator Mode of the On-Board Mis¬ 
sion Planning Function (see Section V of the main report). That is, it deter¬ 
mines how a trajectory will be generated, the subroutines necessary and the 
linking of subroutines. In TGP this logic and prescriptions are dependent on 
the mission type but otherwise fixed. For example, ail missions requiring 
rendezvous with an orbiting target (real or fictitious) are planned in the same 
way; only the results differ. Figure 1-2 is a detailed math flow of MAIN. 

The symbols used may be found in Table 1-1. 

Referring to Figures 1-1 and 1-2, after reading-in and convert¬ 
ing input data MAIN branches if a Surveillance case is to be planned (M-4), 
otherwise the EPHEM sifbrputine is called and certain parameters of the tar- 
igetls-oi bi^vaje computed. If a parking orbit is required (M=3) the LAUNCH 
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routine is used to determine the best time to launch; that is, when the out-of- 

? 

plane angle with the target orbit is minimum. The launch time (tj^) is con¬ 
strained to occur after a reference time (t Q ) and before a maximum time (tj^) 
both of which are inputs. The outputs of LAUNCH are used with UP 1 to 
generate a optional ascent trajectory to a circular parking orbit. 

Once in the parking orbit* the orbital burns necessary to achieve 
rendezvous with the target vehicle are generated as a two-burn maneuver with 
RENDI or as a set of possible three-burn maneuvers with VPG. 

The results from RENDI and YPG are printed-out with the launch 
conditions and ascent phase. Selected options are then ready to be optimized 
with TOP (Appendix II). 

When a direct ascent is called for (M=l, 2), MAIN branches after 
generating the target parameters. The problem here, is to determine the launch 
conditions and an ascent trajectory so that the coast trajectory after cutoff in¬ 
tersects the target orbit when the target is at the intercept point. This is the 
case, for example, when a direct-ascent rendezvous is to be planned. The 
first step is to employ LAUNCH to generate launch conditions for minimum 
out-of-plane angle at the earliest possible launch time (t-^ = t Q ). From this, 
the intercept direction is determined. The EPHEM routine is then called to 
determine the tar get 1 s position and velocity at the intercept point, and the time 
required to reach that point (tf). The intercept position and time to intercept 
(tf - t G ) are inputted to UP 2 which generates an ascent trajectory, if one exists, 
which will arrive at that position at the proper time. If the A V required on 
this trajectory less than the maximum available, a solution has been found 
and the program writes the results and stops. If, on the other hand, the tra¬ 
jectory is not feasible, the program increments the launch time by an input 
constant (Atj^) and proceeds to generate a new trajectory. The launch time is 
incremented in this fashion until a feasible trajectory results (i.e., one with 
A V < AV^^), or a maximum value for t^ ^Lmax' * s exceec ^ ec ^ i* 1 latter 
case, the launch time is fixed at the ,! best ,T time found so far (i.e., corre¬ 
sponding to a minimum value of (AY), and the launch azimuth is now incre¬ 
mented from its minimum value by a fixed amount (A-A^ 


, an inputf. 
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In this iteration the maximum and minimum Values for launch azimuth are set 

I. 

to the desired value in LAUNCH. This results in the launch conditions and 

J 

intercept point when the launch time and the paunch azimuth are specified and 

u 

will,in general, not result in the minimum out-of-plane angle with the target's 
orbit. If a solution can not be found for Aj^< Aj^ max the program stops. 

Within the above mentioned constraints of launch time and launch azimuth, this 
exit indicates that the target cannot be reached with a direct ascent trajectory. 

In the case of surveillance missions (M=4), MAIN initially checks 
the possibility of a direct over-flight by proper cutoff conditions during ascent. 
This is done by determining the maximum inclination attainable from the launch 
site when launch azimuth restrictions are considered. If this maximum incli¬ 
nation exceeds the latitude of the point to be overflown, a first pass overflight 
is calculated. In this application the inputs to UP 1 are varied in an attempt 
to find an elliptic parking orbit which has the desired altitude over the point 
and gives overflight at the desired time. In this iteration the parking orb^t 
parameters and overflight time are varied. 

If a solution cannot be found or if the target 1 s latitude exceeds the 
maximum possible inclination, a parking orbit and transfer maneuver is used. 
UP 1 is used to put the spacecraft into a circular parking orbit which is as 
close to the target as launch azimuth constraints allow. The EPHEM and RVVR 
routines are then used in an iteration loop to determine the time or the orbital 
burn and the overflight time (which defines an aim point for the burn). The 
iterations are continued until a suitable trajectory is found which results in a 
total AV (ascent to parking orbit and orbital maneuver) which is less than the 
vehicle's capability (AV max ). 

1.4.2 LAUNCH SUBROUTINE 

A math flow of the LAUNCH subroutine is presented in Figure 1-3 
and the corresponding symbols are in Table 1-2. This routine is used to de¬ 
termine the launch conditions for: (a) a specified launch time, (b) the launch 
time which results in the minimum possible out-of-plane angle, or (c) a speci¬ 
fied launch time and launen azimuth. 



AI-6 


Referring to Figure 1-3, option (a) corresponds to LCL positive. 

i 

The specified launch time (t) is used to determine the location of the launch 
site relative to the target ! s ascending mode 1^). The launch azimuth (A-^) is 
then computed which will result in the minimum out-of-plane angle with the 
given target orbit plane. This results in a 90 degree range angle between 
launch and target plane intersection (A0 ). The computed launch azimuth is 
compared with the allowed minimum and maximum values an< ^ ^'Lmax^ 

and set equal to the appropriate limit when it is exceeded. If the launch azimuth 
is constrained the range angle A 9 is computed; otherwise it is set to 90 de¬ 
grees. The launch conditions are then specified by generating unit vector 
along the launch radius (rjO, in the launch plane defined by A-^ (v^j and nor¬ 
mal to the launch plane (n^). Finally, the out-of-plane angle between the 
launch and target planes (Ai) is calculated. 

Option (b), when LCL is negative, is used to determine the launch 
times ( Ti and r 2 ) when the minimum possible out-of-plane angle (Ai) can be 
achieved. If the launch site latitude (X^) i s less than the targets inclination 
(i-j,) there are two opportunities each day when A i is zero; otherwise both 
opportunities coincide and result in a value of A i equal to the difference in 
X l and irp. When this option is used, the launch unit vectors r^, v^, n^ are 
obtained by setting t = t \ (or T 2)1 LCL positive, and re-entering LAUNCH. 

The last option (c) is just a special case of (a) where the launch 
azimuth limits (A ! L m j n and A ! j^ max ) are both set equal to the desired launch 
azimuth, A^. 

,1.4.3 UP 1 AND UP 2 SUBROUTINES 

These two subroutines are really just two entry points to one rou¬ 
tine. This is possible because of the large amount of commonality between the 
computations for ascent to a parking orbit (UP 1) and ascent to intercept (UP 2). 
This routine is based on the optimal-explicit guidance algorithm developed 
during this study which is discussed in Par. 6.4 of the main report and de¬ 
scribed in detail in Appendix III. This algorithm is used for the cxoatmo- 
spheric phase of ascent and gives ascent trajectories with nearly minimal Jcmrm 
burning times (and therefore, minimum fuel consumption). 
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In order to initialize the explicit guidance equations, the vehicle T s 
state at the "top of the atmosphere' r (dynamic pressure less than 30 lbs/ft.) 
must be determined in terms of the launch conditions. Presently, both UP 1 
and UP 2 determine these conditions by analytical equations which use three 
vehicle-dependent and payload-dependent parameters. These parameters, 
once determined, are then fixed and used for all launch azimuths. The three 
parameters, in conjunction with the launch conditions, are used to calculate 
the position and velocity at the beginning of explicit guidance by use of spherical 
trigonometry relations. This analytical technique of handling the atmospheric 
phase of ascent is very efficient computationally, but initial simulation re¬ 
sults indicate prediction errors of the order of 1 - 3 seconds in the final in¬ 
jection time. This aspect of planning the ascent phase can and should be im¬ 
proved by a more accurate modeling of the atmospheric phase. 

1.4.4 EPHEM SUBROUTINE 

The ephemeris subroutine EPHEM, see Figure 1-4 and Table 1-3, 
is used to determine the spacecrafts or targets position and velocity (r^ and r) 
at a given epoch based on a given ephemeris (r Q , r Q ) and epoch (T q ). The pro¬ 
gram has two modes of operation: (1) Option A (NOP = 0) when the range angle 
(B) between the epochs is specified, and (2) Option B (NOP = 1) when the trans¬ 
fer time (t) is specified. When Option A is employed, the transfer time (t Q f) 
is calculated in EPHEM in addition to the position and velocity. This routine 
was formulated during the studies of orbit maneuvers and a description of its 
operation and the algorithms employed in it are given in Section VI of the main 
report. 

• / 

1.4.5 RVVR SUBROUTINE 

The Required Velocity Vector Routine (RVVR) is used to determine 
the velocity required to transfer between two positions in space. The slope 
(target of the flight path angle) at arrival at the second point ( m a ) is the con¬ 
straint imposed to completely define the problem. Figure 1-5 is a math flow 
of this subroutine and the symbols are defined in Table 1-4. 
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This routine has five options available for specifying the slope 
and thus, generally, five transfer options (KK = -1, 0, 1, 2, 3). These op¬ 
tions were developed during this study as part of the generalized guidance 
equations subtask in order to study different techniques for generating two- 
burn orbit transfer maneuvers. The rationale and implementation of these 
options is discussed in Section VI of the main report. 

For those cases where a 180 degree transfer is involved and the 
transfer plane becomes undefined, an input quantity, a, is used in the subrou¬ 
tine to define the plane. If this case is not known when the routine is entered, 

0t = 0 and the velocity at the first position is used to define the plane and an 
indicator is set (N180 = +1). The main program can then specify a value for 
a and re-enter the routine, if required. This subroutine calls its own time - 
of-flight routine for calculating the time on the transfer arc: TIME (Figure 
1-6 and Table 1-5) for elliptic arcs and SIMTIM (Figure 1-7 and Table 1-6) 
for parabolic and hyperbolic arcs. 

1.4.6 RENDI SUBROUTINE 

The Rendezvous-Intercept routine (RENDI) is a trajectory planning 
program which is employed to find local minimum - AV two-burns orbit trans¬ 
fers which result in rendezvous with an orbiting target vehicle. The target 
may be fictitious as is the case of placing a payload in a synchronous equatorial 
orbit where a "target 11 is introduced to control the earth longitude of the final 
orbit injection. 

The program uses the EPHEM and RVVR subroutines (paragraphs 
1.4.4 and 1.4.5 above) to generate a two-burn rendezvous and a one-dimen¬ 
sional search routine (Subroutine BEST) to find local minimums of the total AV 
required. In TGP this routine is called to generate all local minimum AV 
rendezvous for which the first burn is allowed to occur between TSTART and 
TSTOP (two input times). 

This routine represents the planning mode of operation of the Ren¬ 
dezvous-Intercept guidance algorithm which is discussed in detail in Section VI 
of the main report. Figuie 1-8ris a math flow of fhls subroutine and Table 1-7 
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contains the symbols employed. The one-dimensional search routine (BEST), 
based on PowelUs technique, is represented in Figure 1-9 with a list of sym¬ 
bols in Table 1-8. 

1.4.7 VPG SUBROUTINE 

The Variable Point Guidance (VPG) subroutine is a planning pro¬ 
gram used to generate a series of orbital burns which result in rendezvous 
with a specified target vehicle. This is an IBM program based on the tech¬ 
nique developed in the Variable Point Guidance and Targeting Study, Refer¬ 
ence ( 1 ). 

This program was developed during this study based on earlier re¬ 
ports generated during the various VPG studies phases (especially references 
1, 8, and 9 of Reference ( 1 ). The program incorporates many of the fea¬ 
tures reported in these earlier papers but differs in actual implementation 
from that reported in Reference ( 1 ). Some of the important features of the 
present version of the VPG subroutine are: 

1. The spacecraft starts from a circular parking orbit. 

2. The target (real or ficticious) is an arbitrary orbit (relative 
to the spacecraft). 

3. All burns are considered as impulses. 

4. Spherical gravity model assumed. 

5. The flight path angle, relative to local horizontal, is un¬ 
changed by a burn (i.e. except for plane change the AV 
added is co-linear with the velocity). 

6. All burns except the first occur on the iine-of-nodes between 
the spacecraft and target orbits. 

7. The first burn is positioned in the parking orbit to satisfy 
(5) and (6). 

8. Eight options are generated for each input case. Four op¬ 
tions correspond to rendezvous occurring at one node, the 
other four to rendezvous at the opposite node. The four op- 

- 5 d.-V-tion-s are Bielliptic Chase, Bielliptic Lob, Full Orbit Phasing 
Chase, and Full Orbit Phasing Lob. 
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9. Bielliptic maneuvers involve three major burns: the first 
at parking orbit altitude, the last at target altitude (on the 
line-of-nodes) and the second 180 degrees from the third. 

The altitude of the second burn is adjusted to insure proper 
phasing at the third burn. If this altitude is below the target 
orbit, it is a Chase solution, if above, it is a Lob solution. 

10. Full Orbit Phasing maneuvers also involve three major burns: 
The first at the parking orbit and the second and third coin¬ 
cident at the target altitude. The second burn, at the target 
orbit, is adjusted to produce a phasing orbit such that after 
an integer number of revolutions the spacecraft and target 
are coincident. The third burn then matches the spacecrafts 
velocity to that of the target. If the period of the phasing 
orbit is less than that of the target, it is a Chase solution, 
otherwise it is a Lob solution. 

11. The plane change to be effected at each burn is determined 
by optimally splitting the total change between the major 
burns. This routine includes an option (input) of also con¬ 
sidering a pure plane change maneuver at target altitude. 

One of the major differences between this program and that in 
Reference ( 1 ) is that the ascent to orbit is not included in the VPG subrou¬ 
tine. This phase of the mission is generated with the UP 1 subroutine ( Para¬ 
graph 1.4. 3 of this Appendix) because of the flexibility and performance it offers. 

A detailed description of the philccophy underlying the Variable 
Point Guidance and Targeting technique as well as its implementation can be 
found in Reference ( 1 ). A math flow of the VPG subroutine can be found in 
Figure I-10 and the symbols are listed in Table I-9. VPG itself employs 
three subroutines: an optimal Plane Change Angle routine - PCA (Figure I-11 
and Table I-10); a Required Velocity routine - REQ (Figure 1-12 and Table 
I—11); and a Time-of-Flight routine - TOF (Figure 1-13 and Table 1-12). 
The latter two .routines (REQ and TOF) were generated when VPG was initially 
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altitude (LCO = +1, LCS = +1) the radius magnitude is also calculated inter¬ 
nally to insure a 24 hour period (R^ is then used only to define the target ! s 
direction at*the epoch, t*p). 

If the mission is to start from a parking orbit (LCL = -1) the space¬ 
crafts ephemeris is inputted: 

Rp = Parking orbit position at t = t 

Vp = Parking orbit velocity at t = t 

t = Epoch time. 

P 

The type of mission to be planned is determined by the integer M: 

M = 1 Intercept 

M = 2 Direct Ascent Rendezvous 

M = 3 Orbital Rendezvous 

M = 4 Surveillance - from parking orbit 

M = 5 Surveillance by first pass overflight (set in program if 

found possible) 

For M = 4, the surveillance point is defined by: 

\ T = Latitude 

fl ^ = Longitude relative to launch site 

R^ = Overflight radius 

= Error tolerance allowed in the overpass. 

The additional input data required are: the circular parking orbit 
radius Rp (when LCL = +1); last possible launch time last possible 

final time tf max ; maximum AV capability of the vehicle an d incre¬ 

ments of launch time and launch azimuth (Atj^ and AAj^) which are used when 
M = 1 or 2 to find a feasible direct ascent trajectory. 

1*6 TEST CASES 

Three test cases are presented in order to illustrate the perform¬ 
ance of TGP* Table T >13 contains a description of the cases, and Table 1-14 
illustrates the results obtained with TGP. 
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written as an independent program. They were retained when VPG was con¬ 
certed into a subroutine for TGP (instead of .using RVVR and EPHEM) because 
of the internal structure of VPG which would have required a major overhaul 
in order to make use of the more general RVVR and EPHEM routines. 

1 • 5 PROBLEM SPECIFICATION 

The coordinate system used to define the problem and in which all 
calculations are performed is the XYZ ECI (Earth Centered Inertial) system, 
illustrated in Figure 1-14. The XY plane is the equatorial plane and the Z- 
axis is through the north pole. This figure also illustrates some of the param¬ 
eters involved in the program, which are defined in Table I-l. 

If the mission is to be planned from launch (LCL = +1) the launch 
site is specified by: 

X = Launch site latitude. 

L 

,0, _ = Launch site longitude relative to X-axis at t 

Lo o 

t = Reference time, 

o 

A* . = Minimum allowable launch azimuth. 

Lmm 

A* = Maximum allowable launch azimuth. 

Lmax 

Two alternative ways are provided for specifying the target orbit 
If LCO = -1, the target orbit is specified by its ephemeris at some epoch: 

R = Targets position at t = t^ 

V = Targets velocity at t = t^ 

= Epoch time; or, 

If a circular target orbit is desired (LCO = +1) the orbit is specified by: 

n = Unit normal to target orbit 

c 

R = Targets position at t = t^ 

t^ = Epoch time. 

and the program computes the circular velocity from !and'form3,the ver vr 
locity vector with n^ and R»p. When the target orbit corresponds to a synchronous 
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Case 1 - This case requires planning the launch conditions and 
orbital maneuvers to achieve rendezvous with a fictitious target in a synchro- 
nous equatorial orbit stationary at launch site longitude. The program gen¬ 
erates, with the LAUNCH subroutine, a launch at tj^ “ 0 at an azimuth of 90° 
which gives the minimum out-of-plane angle possible (A* = 28.5617° = Xj^). 

The ascent to the circular parking orbit is calculated with UP 1 and requires 
483.7 seconds with 14, 444 feet/sec. of A V left for orbit maneuvers. 

The orbit transfer necessary for rendezvous is generated as a two- 
burn maneuver with RENDI and as a series of three-burn maneuvers with VPG. 

In the case of RENDI, Table 1-14 lists the results of the first local minimum 
A V found and the best local minimum in the time span allowed for the search 
(t£ < 1 day). The VPG results in Table 1-14 represent one Bielliptic and one 
Full Orbit Phasing Option. Both are three-burn maneuvers and represent 
maneuvers which require less Av than is available (14,444 feet/sec.) and 
which, if possible, achieve rendezvous within 1 day (86,400 seconds). 

The Bielliptic option requires 14, 138 feet per second but final 
rendezvous is at 102, 716 seconds. This is accomplished by gross phasing in 
the parking orbit for 7 full orbits after the first nodal crossing. The rendez¬ 
vous is then effected by two 180 degree transfers, with the second burn being 
above the target orbit at 143, 437, 670 feet (target radius is 138, 607, 380 feet). 

The second VPG option is a Full Orbit Phasing case which requires 
14,032 feet/second of AV with rendezvous at 64, 800 seconds. In this case, the 
phasing is accomplished with one orbit at initial altitude followed by a 180 de¬ 
gree transfer to target altitude. The second Durn then p^ts the spacecraft into 
a phasing ellipse with a period of 38,990 seconds (target period is 86,400 sec¬ 
onds) and after one phasing orbit the third burn matches the velocities for ren¬ 
dezvous. 

Case 2 - The target in this case is at Space Station Altitude (185 
nm) and the problem is to plan a rendezvous mission from launch. The LAUNCH 
routine is used to find the most favorable launch conditions, which are t^ = 
i44* 502 seconds, and Aj^ = 74. 9178 degrees. This results in the ascent trajec¬ 
tory intersecting the target orbit 90 degrees down range with zero out-of-plane 
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angle (irj. > X so an in-orbit launch time opportunity is available twice per 
day). The ascent to a 100 nm parking orbit is then generated with UP 1 and 
results in an injection time of 15, 089 seconds (585 seconds after launch) with 
8, 762 feet/sec. of AY remaining for the rendezvous maneuvers. 

The orbital maneuver represented in Table I -14 is a VPG Full 
Orbit Phasing option which requires 6, 057 feet/sec. of AV with rendezvous 
at 27, 383 seconds (12, 881 seconds after launch). The maneuver involves a 
180 degree transfer to target altitude followed by one revolution in a phasing 
orbit. The period of the phasing orbit is greater than tha.t of the target orbit 
(8,280 seconds vs. 5,485 seconds) so that this is a Lob option. 

Case 3 - The final test case is the same as Case 2 except that a 
direct ascent rendezvous is to be planned. TOP finds a solution corresponding 
to the launch site 80. 25 degrees beyond the target’s ascending node (tj^ = 19,260 
seconds);this corresponds to approximately 80 minutes after the launch site 
has passed through the target orbit plane. The launch azimuth of 84.8443 de¬ 
grees results in the intercept point being 90 degrees down range and an out- 
of-plane angle of 3.05 degrees. These conditions are inputted to UP 2 with 
the time that the target will arrive at the intercept point (2112 seconds after 
launch). An ascent trajectory is generated which achieves the desired final 
state and has 2, 734 feet/sec. of AV left after the adaptation burn for rendez¬ 


vous. 
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TABLE I-lj 

MAIN PROGRAM SYMBOLS 


Math Flow 
Symbol 

Description 

Program 

Symbol 

M 

Run Type Variable 

M 

t. 
o 

Reference Time 

TREE 

Lo 

Initial Position Vector (Target) 

RTO 

Zto 1 

Initial Velocity Vector (Target) 

VTO 

l TO 

Initial Time (Target) 

TTO 

/“ 

Gravitational Constant = 1.4082878x10^ 

GC 

r 

Position Vector 

RTN 

V 

Velocity Vector 

VTN 

TOE 

Time of Flight 

TOF 

MA 

Slope at Arrival 

SLT 

a T 

Semi Major Axis of Target Orbit 

1 

A 

h T 

Angular Momentum 

ANG 

-St 

Unit Cross Product Vector (Target) 

UNT 

h T 

Target Semi Latus Rectum 

SEMI.UR 

®T 

Eccentricity 

E 

X T 

Target Inclination Angle 

EYE 


Longitude of Ascending Node 

CLONG 

al max 

Maximum Allowable Launch Azimuth Angle 

ALMAX 


(Sheet 1 of 5) 
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Table 1-1 Main Program Symbols (Continued) 


Math Flow 
^ Symbol 

Description 

Program 

Symbol 

al min 

Minimum Allowable Launch Azimuth Angle 

ALMIN 

T 1 

1st Zero Out-of-Plane Launch Time 

TAU1 

T 2 

2nd Zero Out-of-Plane Launch Time 

TAU2 

X L 

Launch Site Latitude 

SLAT 

12 LO 

Launch Site Longitude at t 

SLONG 

t 

Desired Launch Time 

T 


Launch Azimuth 

AL 

A i 

Out-of-Plane Angle 

DELEYE 

A0 

Range Angle to Target Plane Intersection 

TANG 

-i.L 

Unit Vector Thru Launch Site 

RL 

C _r L 

Unit Vector In Launch Plane 

VL 

-2l 

Unit Vector Normal to Launch Plane 

UNITL 

N 

Program Counter = Number of Launch Time 



Calculations 

N 

t l, 2, 3,n 

Launch Times 

TP(N) 

*1 

Best Launch Time 

T i 

*1 

Unit Vector Thru Launch Site at t^ 

R i 

r e 

Earth Radius 

RE 

_?1 

Unit Vector Normal to Launch Plane at t^ 

UNIT 1 

r 

a 

Apogee Radius 

RA 

r 

P 

Pergee Radius 

RP 

f o 

Parking Orbit Altitude -L 

' RPARK 


Launch Site Longitude at Launch Time 

XLONG 


(Sheet 2 of 5) 






AI-18 


Table 1-1 Main Program Symbols (Continued) 


Math Flow 
Symbol 

Description, 

• Program 
Symbol 

AV* 

Current Max Value of AV-^ 

D VS TAR 

CD 

Range Angle to Target 

TH.ETF 

V 

Time of Insertion or Rendezvous 

TFP 

3, 

Position Vector at t^ T 

RF 

Jf 

Velocity Vector at t ^ 

VF 

t * 

Time of Conditions with Highest Remaining AV 

TSTAR 

a l* 

Launch Azimuth with Highest Remaining AV 

ALSTAR 

LC1 

Program Indicator to Step Launch Azimuth 

PCI 

aa t 

Launch Azimuth Increment 

AZ DEL 

.. At i 

Launch Time Increment 

TDEL 

t l MAX 

Latest Possible Launch Time 

TLMAX 

*MAX 

Maximum Inclination Angle 

EYEMAX 

XT 

i 

Target Latitude 

TLAT 

AT 

Elapsed Time 

DELT 

*4 

Transfer Time for Zero Flight Path Angle 

T4 


Overpass Position Vector 

RT 


Time at Overpass Position 

TF 

6 f 

Error in Time at Overpass Position 

ERRF 

4 fMAX 

Last Possible Final Time 

TFMAX 
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Table 1-1 Main Program Symbols (Continued) AI-19 


^ Math Flow 

^ Symbol 

Description 

9 

Program 

Symbol 

W E 

Earth Rotational Rate 

WE 

‘2 

Time at Injection into Parking Orbit 

T2 

*2 

Parking Orbit Position or Cut-Off Vector at Time t 

Ca 

R2 

Iz 

Parking Orbit or Cut-Off Velocity Vector at Time t 

V2 

AV 

M 

Remaining AV(Ft/Sec) Available 

DELVM 

av r 

Remaining AV(Ft/Sec) Available 

DEL.VR 

Ji 

Target Position Vector at 

RTN 

Jff 

Target Velocity Vector at t 

VTN 

n- 

Unit Cross Product Vector of R_ and V 

CP 

— 2 

—•’2 —2 


^ a 

Cosine of Out-of-Plane Angle Between Parking 



Orbit and Target Orbit 

ALPHA 

n_ T 

i 

Unit Vector to Line-of-Nodes 1 

CPPR 




e 

o 

Transfer Angle from Parking Orbit to Rendezvous 

i 

THE TO 

*3 

Position Vector at Rendezvous 

R3 

J 3 

Velocity Vector at Rendezvous 

V3 

s 

1 

Time of Rendezvous 

i 

T3 

— T3 

Target Position Vector at Rendezvous 

RT3 

V 

Target Velocity Vector at Rendezvous 

VT3 

—T3 



A£ 3 

Angular Position of Target Relative to Line-of-Nodes 



at t 3 

DELF3 

f 3 

True Anomaly of Target at t 

J \ 

F3 

r £ ln ; -'" y 

‘ Angular Position of Line-of-Node Relative to Targets 


Lv. 

Perigee 

FLN 
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Table I-1 Main Program Symbols (Continued) 


Math Flow 
Symbol 


Description 


Program 

Symbol 


Velocity Required for Overflight 


DELTV3 


Angular Overpass Error 


DEL 


Maximum Allowable Overpass Error 


DELSTR 


Final Launch Longitude 


XLONGF 


Final Target Longitude 


TLONGF 


Target Longitude Relative to Launch Site 


TLONG 


Unit Vector to Parking Orbits Descending Node 


CPDN 


Remaining AV 


DELVN 


Range Angle Error 


DELONG 


LCO 


Parameter to Select Circular Orbit Option 


LCO 


Parameter to Select Synchronous Orbit Option 


LCS 


LCL 


Parameter to Select Parking Orbit Data 


LCL 


LCT 


Parameter to Select VPG Option 


R 

— po 

V 

— P° 


Parking Orbit Position Vector 
Parking Orbit Velocity Vector 
Parking Orbit Epoch Time 


RPO 


VPO 


TPO 


Radius of Circular Orbit 


Velocity of Circular Orbit 
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TABLE. ■ 1-2 


AI-21 


LAUNCH SUBROUTINE SYMBOL TABLE 


Math 

Symbol 

© 

Description 

FORTRAN 

Program 

Symbol 

LCL 

Option Indicator: +1 when launch time is specified; 

-1 determines best launch time(s) 

LC 

X L 

Launch site latitude (degrees) 

LA ML 

^LO 

Launch site longitude relative to x-axis at t^ 
(degrees) 

OLO 


Longitude of ascending node of target orbit 

(degrees) 

OT 

l T 

Inclination of target orbit plane to x-y plane 
(i^ > 0) (degrees) 

IT 

t 

Desired launch time 

T 

t o 

Reference time (sec) 

TO 

A ! 

L max 

Maximum permissable launch azimuth (degrees) 

ALMAX 

A ! t • 

L mm 

Minimum permissable launch azimuth (degrees) 

ALMIN 

T i 

Time at which launch site is nearest target plane 

T 1 

T 2 

Time at which launch site is nearest target plane 

T 2 

a l 

Launch azimuth (degrees) 

AL 

Ai 

Out of plane angle between launch plane and target 

orbit (degrees) 

DELI 

a e 

Range angle from launch site to intersection with 

target orbit 

DELTH 

.Jl 

Unit vector through launch site 

RL 


c 
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Table 1-2 Launch Subroutine Symbol Table (Continued) 


Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

V 

— L 

Unit vector in the launch plane and perpendicular 

to r 

RL 

n 

—. L 

Unit vector normal to launch plane 1 

NL 


< i.L X J L > 


(x) 

E 

Rate of the earth’s rotation (rad/sec) * 

WE 


Longitude of launch site at time t relative to 

target orbit’s ascending node 

PHIL 

*IJ 

Value of (p when t = r ^ 

PHIL1 

^L2 

Value of <f> when t = TV 

PHIL2 

LL 2 


„ A ia 

Value of A , at t = T , for minimum Ai 

L 1 

AL1 

V 

> 

It* 

CO 

Value of A , at t = T , for minimum Ai 

Li 2 

AL2 

Ai i 

Value of Ai at t = 

DELI1 

Ai 

Value of Ai at t = T_ 

DELI2 

2 

2 
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TABLE 1-3 ’ 

EPHEM (EPHEMERIS SUBROUTINE) 
SYMBOL TABLE 


• 


FORTRAN 



Program 

Math ; 

; 

Description 

Symbol 

Symbol 


, 

r , r* 

_o __o 

Initial position and velocity (ft. and ft/sec) 

RO, RODOT 

TO 

Initial time (sec) 

TO 

P- 

Gravitational constant 

U 

e 

Transfer angle (radian) 

TH 

t 

Time of transfer-input for option B (sec) 

T 

DLE 

Convergence criteria 

DLE 

NOP 

= 0 for Option A (angle input) 

= 1 for Option B (time input) 

NOP 

r , r 

Position and velocity at end of transfer 

R, RDOT 

tof 

Time of flight (Option A only) (sec) 

TOF 

m 

a 

Arrival slope 

MA 

a 

Semi-major axis (ft) 

A 

r 

o 

Magnitude of 

ROM 

e ro 

Unit radial vector 

ERO 

V 

ro 

Initial radial velocity (ft/sec) 

VRO 

C 

2 

Moment of momentum (ft /sec) 

C 

C 

2 , 

Magnitude of C (ft /sec) 

CM 

V So 

Initial horizontal velocity (ft/sec) 

VTHO 

m 

o 

Initial slope 

XMO 

v r 

Magnitude of r - 

RA 

AE 

Eccentric anomaly (radian) 

DLEN 
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Table 1-3 EPHEM (Ephemeris Subroutine) Symbol Table (Continued) AI-24 




FORTRAN 

Math 


Program 

Symbol 

• Description 

Symbol 

p 

Period (sec) 

p 

V 

ra 

Radial velocity at arrival 

VRA 

V 0a 

Horizontal velocity at arrival 

VTHA 
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TABLE 1-4 


AI-25 


RVVR (REQUIRED VELOCITY VECTOR SUBROUTINE) 


SYMBOL TABLE 


Math 

Symbol 

Description 

FORTRAN 

1 

Program 

Symbol 

u 

o 

o ‘ 

Initial position and velocity (ft. and ft/sec) 

RO, RODOT 

r 

a 

Aim point (ft) 

RA 

m 
. a 

Slope at arrival 

XMA 

X kl 

o 

Constant to determine 180 transfers 

XK1 

P- 

Gravitation constant 

U 

NN 

Number of complete orbits during transfer (input) 

NN 

r Ta 

Magnitude of target radius (Option A-3 only) 

RTAM 

KK 

= -1 for Option B of Rendezvous -Intercept 

KK 


0 for Option A -1 of Rendezvous-Intercept 



1 for Option A-2 of Rendezvous-Intercept 



2 for Option A-3 of Rendezvous-Intercept 



3 for Option A-4 of Rendezvous-Intercept 


ttl 

Po 

Slope at r^ before velocity change 

XPO 

V 

Unit horizontal vector 

ETHR 

a 

Plane change angle used in 180° transfers (radian) 

ALF 

N180 

= 1 when 180° transfer has been performed 

N180 


= 0 otherwise 


e 

a 

Transfer angle (radian) 

THA 

tof 

Time of flight (sec) 

TA 

V 

-» r 

Velocity required for transfer (ft/sec) 

VR 

NSIG 

Indicates type of conic transfer 

NSIG 

a 

Semi-major axis (ft) > ior 

: . . . . . "j 

A, AX 
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Table 1-4 RYYR (Required Velocity Vector Subroutine) Symbol Table (Continued) AI-26 


Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

r 

Magnitude of r (ft) 

ROM 

o 

— ► o 


r 

Magnitude of r (ft) 

RAM 

a 

— a 


e 

Unit radial vector 

ER 

r 



r 

Radial component of r (ft) 

RAR 

— ar 

a 


Ja6 

Horizontal component of r^ (ft) 

RATH 

r a0 

Magnitude of r q (ft) 

RA.THM 

\< 

CD 

Horizontal component of r (ft/sec) 

o 

VTH 

v e 

Magnitude of Vg (ft/sec) 

•* _ i 

VTHM 


c 


c 


(Sheet 2 of 2) 



TABLE 1-5 


TIME SUBROUTINE SYMBOL TABLE 




Math 


Symbol 


Description 



cos 0 a 
sin 6 a 
m 

a 

r 

a 

NN 

tof 

C 

a 

AE 


Horizontal component of velocity (ft/sec) 
Radial component of velocity (ft./sec) 
Gravitational constant 
Magnitude of position vector (ft) 

Transfer angle (radian) 


Slope at arrival 

Magnitude of arrival point radius (ft) 
Number of complete orbits 
Time of flight (sec) 

Magnitude of angular momentum 
Semi-major axis (ft) 

Eccentric anomaly (radian) 


AI-27 

FORTRAN 

Program 

Symbol 

VTHR 

VRR 

U 

ROM 

CTHA 

STHA 

XMA 

RAM 

NN 

TA 

CA 

A 

DE 




TABLE 1-6 


AI-28 


SIMTIM SUBROUTINE SYMBOL TABLE 


Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

e 

Transfer angle (radian) 

TH 

m 

o 

Departure slope 

XMO 

r 

p 

Departure radius (ft) 

ROM 

C 

Magnitude of angular momentum 

H 

m 

a 

Arrival slope 

MA 


Gravitational constant 

U 

tof 

! 

Time of flight (sec) 

TOF, TOFN 


( 


c 


TABLE 1-7 


AI-29 


( 

Math . 
Symbol 

RENDI - RENDEZVOUS-INTERCEPT SUBROUTINE 

SYMBOL TABLE 

Description 

FORTRAN 

Program 

Symbol 

r 

—• i 

Initial interceptor position (ft) 

RI 

• 

r. 

1 

Initial velocity of interceptor (ft/sec) 

RIDOT 

-iTi 

Initial target position (ft) 

RTI 

-iTi 

Initial velocity of target (ft/sec) 

RTIDOT 

TSTOP 

Maximum time at which first burn can occur (sec) 

TSTOP 

m 

a 

Arrival slope 

MA 

k i 

c 

= -1 for RENDI Option B 

= 0 for RENDI Option A-l 

= 1 for RENDI Option A-2 

= 2 for RENDI Option A-3 

= 3 for RENDI Option A-4 

Kl 

NOIT 

= 0 return when t + TSTART > TSTOP 

= 1 return after first minimum AV is attained or 

t + TSTART > TSTOP 

= -1 return when AV < AV , or 

M 

t + TSTART > TSTOP 

NOIT 

TSTART 

• • 

Reference time when r., r., r m ., and r m . are 

— l — l — Ti —Ti 

given (sec) 

TSTART 

AV 

M 

Maximum allowable AV for trajectory (ft/sec) 

DELVM 

H- 

Gravitational constant 

U 

AT 

Amount by which the time of first burn is incremented 

in search for minimum AV (sec) -o 

DELT 
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Table 1-7 RENDI - Rendezvous-Intercept Subroutine Symbol Table (Continued) 


Math 

Symbol 


Description 


FORTRAN 

Program 

Symbol 


r , r 
o o 


Time of first burn relative to TSTART (sec) 

Interceptor position and velocity at time of first 
burn 


RO, RODOT 


r . r m 

— T—T 


Target position and velocity at time of first burn 


RT, RTDOT 


r. 

inp 


DLIM 


Projection of r^ onto target plane (ft) 

Angle between r. and r m (radian) 

° —'inp — T 

Lower limit for transfer angle (radian) 


RINP 


THH 

DLIM 


ULIM 


Upper limit for transfer angle (radian) 


ULIM 


Period of interceptor (sec) 

Past guess at transfer angle (radian) 
Present guess at transfer angle (radian) 


r TA’ r TA 


Target position and velocity after being transferred in RTA, RTADOT 


orbit through angle 9^ (or 9^) from r, 


Target slope at r r 


MTA 


T or 
TAO 


Time of target transfer from r to r resulting 
from transfer angle of 9 or 9 (sec) 


TTAO or TTA1 


Aim point (ft) 

Desired interceptor arrival slope 


T A or 
AO 


Time of interceptor transfer from r to r correS' 

-»o _a 

ponding to target transfer angle 9 ^ or 9 (sec) 


TAO or TA1 


VALUE 


The function which is reduced to zero by regula 
falsi iterations (sec) 


VALUE 
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Table 1-7 RENDI - Rendezvous-Intercept Subroutine Symbol Table (Continued) 


Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

DERIY 

Approximated derivative of VALUE 

DERIY 

e 2 

New value of Q (radian) 

TH2 

V 

_► r 

Velocity required by interceptor a.t r to transfer 

—► o 

to r (ft/sec) 

a 

VR 

t a 

Fixed transfer time (sec) 

TA 

m A0 

Past guess at arrival slope 

MAO 

m Al 

Present guess at arrival slope 

MAI 

m A2 

New guess at arrival slope 

MA2 

V 

— a 

Arrival velocity of interceptor (ft/sec) 

RADCT 

L AV 

Total change in velocity (ft/sec) 

DELV 
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TABLE 1-8 


BEST-SUBROUTINE SYMBOL TABLE 


^ Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

DY 

Velocity change (AV) 

DV 

T 

Time of first burn 

T 

TH 

Transfer angle 

TH 

IPL 

Counter to tell when to apply quadratic fit 

IPL 

DT 

Time increment equal to 1/4 the interceptor or 

target period, whichever is smaller 

DT 

T MINS A 

T at last minimum found 

T MINS A 

NOPT 

= 1 when minimum has been found 

NOPT 


= 0 otherwise 


£ 8 T 

Time increment = DT/10 

SDT 

V 

TM(I) 

1=1,4 

Storage for first burn times 

TN(I) 

DVEL(I) 

1=1,4 

Storage for AV ! s corresponding to TM(I) 

DVEL(I) 

1=1, 4 

Normalized values of TM(I) 

TM(I) 

T 

Normalized value of time indicated by quadratic to 



give minimum AV 




TABLE 1-9 


AI-33 



VPG - VARIABLE POINT GUIDANCE SUBROUTINE 



SYMBOL TABLE 


Math 

Symbol 

* 

Description 

FORTRAN 

Program 

Symbol 


I. Inputs 


a T 

Target semi-major axis (ft) 

AT 

®T 

Target eccentricity 

ET 

f LN 

Angular Separation of Line-of-Nodes and Target's 

Line-of-Apsides (deg) 

FLN 

Af 3 

Position of target relative to Line-of-Nodes at t 

(deg) 

DELF3 

t 

3 

Time of spacecraft’s first nodal crossing (sec) 

T3 

R 

P 

Radius of spacecraft’s circular parking orbit (ft) 

RP 

Ai 

Out-of-plane angle between spacecraft (S/C) and 

target orbits 

DELI 

n 

— Z 

Unit vector normal to S/C parking orbit 

EN2 

_J?N 

Unit normal along Line-of-Nodes (defines Node A, 

ENN 

Node B defined by - n ) 


H- 

i 

Gravitation constant 

U 


II. Internal Symbols 


AR 

max 

Maximum altitude above target orbit allowed for 

intermediate aim point in Bielliptic option (ft) 

DELRM 

R* . 

mm 

Minimum radius allowed for intermediate aim point 

in Bielliptic option (ft) 

RMINS 

AR *■ 
o 

Initial altitude increment, added to parking orbit 

DELROS 

radius, used for intermediate aim point in Bielliptic 



" option (ft) 
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Table 1-9 VPG - Variable Point Guidance Subroutine Symbol Table 
(Continued) 


AI-34 


Math ;> 

Symbol 

Description 

FORTRAN 

Program 

Symbols 


Error tolerance used to stop iterations in optimum 

EPA 


plane change angle routine (rad) 


S Ai 

Sign of the Ai input to subroutine 

SDEL1 

"P 

Angular rate of S/C in its orbit (rad/sec) 

WP 

P P 

Period of S/C orbit (sec) 

PP 

V P 

Velocity in circular parking orbit (ft/sec) 

VP 

P T 

Period of target orbit (sec) 

PT 

h T 

2 

Angular momentum of target orbit (ft /sec) 

HT 

r ta 

Target radius at Node A (ft), see n^ 

RTA 

r tb 

Target radius at Node B (ft), seejn^ 

RTB 

R T3 

Target radius at t = t^ 

RT3 

V TA 

Target ! s velocity at Node A (ft/sec) 

VTA 

V TB 

Targets velocity at Node B (ft/sec) 

VTB 

Y 

T3 

Target 1 s velocity at t = t 

VT3 

y ta 

Target’s flight path angle at Node A (rad) 

GAMTA 

y tb 

Target’s flight path angle at Node B (rad) 

GAMTB 

R. i = 3, 7 

1 

S/C radius at t = t., i = 3, • - • , 7 

1 

R3, • • • , R' 

V.~ i = 3, 7 

S/C velocity at t = t. , i = 3, • • • , 7 

V3M, 

l 

1 

and VM(I) I 

V. + , i = 3, 7 

| 

S/C velocity at t = t. , i = 3, • • • , 7 

V3PL, 

1 

1 

1 . 

and VPL(I) 

R* 

^ max 

Makimum radiifb'allowed foP intermediate, aim. point 

RMAXS 


in Bielliptic option 

■ 
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Table 1-9 

^ Math : 

Symbol 

1 VPG - Variable Point Guidance Subroutine Symbol Table 

(Continued) 

* Description 

AI-35 

FORTRAN 

Program 

Symbol 

A T 

Time of flight (sec) 

DELT 

T 

N 

Time for target to reach rendezvous point for the 

TN 


first time (sec) 


R* 

Radius of intermediate aim point in Bielliptic option 

RS 


(ft) 


m* 

Slope of S/C orbit at R* 

MS 

m 6 

Slope of target orbit at rendezvous point 

M6 

V 

r 

Radial component of required velocity (ft/sec) 

YR 

v <? 

Tangential component of required velocity (ft/sec) 

VTH 

X,i = 3,7 

Fliight path angle of S/C orbit at t = t., i = 3, 7 j 

GAM1, . . . , and 

f 

1 

j 

GAM(I) 1 = 3, 7 

•; t 

T1 

! 

Transfer time to intermediate aim point (sec) 

TT1 

| A0 

■ 

Position of first orbital burn (t = t^) relative to iine- 

DELTH 


of-nodes (rad) 



Flight path angle at R* 

gams 

T 

T2 

Transfer time on second arc in Bielliptic option (sec) 

TT2 

T_ 

T 

Total transfer time, relative to ty for S/C excluding 

TT 


integral number of phasing orbits 


X' 

Non-integer number of full parking orbits (for phasing) 

XP 


before first burn 


A R* 

Increment to R* to correct phasing error (ft) 

DELRS 


Initial value of P__ (sec) 

PPEO 

PEO 

PE 


; P PE 

1 C 

Period of phasing ellipse in Full Orbit Phasing option 

PPE 

4 ' 

I 

(Sheet 3 of 4) 




Table 1-9 V 

(« 

s ' Math 

Symbol 

PG - Variable Point Guidance Subroutine Symbol Table 

Continued) 

Description 

AI-36 

FORTRAN 

Program 

Symbol 

a 

PE 

Semi-major axis of phasing ellipse 

APE 

a. i = 3, 7 

l 

Plane change angle at t = t. i = 3, • • • , 7 

ALFA3, ••• and 

i 

AF(I) I = 3,7 

AY. i = 3, 7 

Velocity impulse at t = t i = 3, • • • , 7 

DV3, • • • , DV7 

and DELV(I) 1=3,7 

n. i = 2, 7 

— l 

Unit normal to plane defined by R. and V. 

— l —l 

EN( ,I),I = 2,7 

e. i = 3, 6 

Unit vector along position vector at t = h i = 3,• • • , 6 

EVEC ( , I), 1=3.6 

R i = 3, 6 

Position vector at t. i = 3, ••• , 6 

l 

RVEC ( , I), 1=3, 6 

_3vi ‘ = 3 - 6 

Unit vector normal to R. and n. 

— l —l 

EVC ( ,1), 1=3, 6 

_V., i = 3, 6 

Velocity vector at R.. 

VVEC ( , I), 1=3,6 

AY., i = 3, 6 

— i 

Total velocity impulse at t. i = 3, • • • , 6 

DLVL ( ,1), 

I = 3, 6 

I 

Integer, full orbits in parking orbit 

I 

J 

Integer full orbits in phasing orbit (Full Orbit Phasing 

option) 

J 

K 

Integer, full orbits of target before rendezvous 

K 
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TABLE 1-10 


PCA - PLANE CHANGE ANGLE SUBROUTINE 
SYMBOL TABLE 


AI-37 


r 

^ Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

J, K, L 

Controls for indices on various vectors 

J, K, L 

NM 

Iteration Counter 

NM 

Ai 

Total plane change angle 

DELI 

v'V) 

Velocity magnitude before I*k burn 

YM (I) 

V + (I) 

Velocity magnitude after l^ a burn 

VPL (I) 

r a) 

Flight path angle during burn 

GAM (I) 

a (I) 

Plane change angle at 1^ burn 

ALFA (I) 

AV (I) 

O 

Impulse at 1^ burn with no plane change 

DELVO (I) 

AY (I) 

Total impulse at I& 1 burn 

DELV (I) 




TABLE 1-11 


AI-38 

REQ SUBROUTINE 
SYMBOL TABLE 


Math 

Symbol 

Description 

FORTRAN 

Program 

Symbol 

r i 

Radius of vehicle position (ft) 

R1 

r 2 

Radius of aim point (ft) 

R2 

e 

Angle between vectors to vehicle position and aim 

TH 


point (radian) 


m 

Slope at arrival at aim point 

M 

H- 

Gravitational constant 

U 

v a 

Horizontal component of required velocity (ft/sec) 

VTH 

Y 

r 

Radial component of required velocity (ft/sec) 

VR 


c 





TABLE 1-12 


AI-39 


TOF SUBROUTINE 
SYMBOL TABLE 


Math 

Symbol 

Description 

FORTRAN 
/ Program 
Symbol 


Radius of vehicle position (ft) 

R1 

r 2 

Radius of aim point (ft) 

R2 

a e 

Angle between vectors to vehicle position and aim 

point (radian) 

DELTH 

v i 

Velocity of vehicle (ft/sec) 

VI 

H- 

Gravitational constant 

j 

U 

r, 

Slope of vehicle 

GAMl 

AT 

Time of flight 

DELT 

! 

h 

Magnitude of angular momentum 

■ 

H 

a 

Semi-major axis (ft) 

A 

m 

Slope at arrival at aim point 

XM 

AE 

Eccentric anomaly (radian) 

DE 


c 



Table 1-13 


TEST CASES 


C9-S6 

Initial 

Conditions 

Final 

Orbit 

Comments 

i 

Launch Latitude: = 28. 5617° 

Semi-Major Axis: 

Rendezvous with fictitious .target 


Launch Longitude: XI = 0 

Reference Time: t Q = 0 

Launch Azimuth Limits: 65°<A T < 
115° L 

a T = 138, 607, 380 ft. 

in synchronous, equatorial 


Eccentricity: e^ = 0 
Inclination: ir^ = 0 

orbit at launch site longitude. 

2 

Same as Case 1 

Semi-Major Axis: 

a T = 22, 060, 185 ft. 
Eccentricity: ~ 0 

Inclination: irp = 32° 

Longitude of Ascending Node: 

Rendezvous with vehicle in 185 
n. rn. circular parking orbit. 



o 

!i 

H 

c5 


3 

Same as Case 1 but no parking 

orbit used, and XI - 80° 

LO 

Same as Case 2 

■ 

- 

Direction Ascent Rendezv6us to 
target of Case 2 


< 


AI-40 

















Table 1-14 

RESULTS QF TEST CASES 


Case 

Launch 

Ascent 



Orbita 

1 Maneuvers 




® Launch Time 

CD 

o 

rtf 

s 

•H 

N 

< 

rtf 

o 

tf 

tf 

d 

(deg) 

r — 

0- Out-Of-Plane Angle 

2. With Target Orbit 

co Cutoff Time 

CD 

O 

9" Altitude At Cutoff 

• 

g (Parking Orbit for 
.L, Cases 1 and 2) 

{$ AV Available for 
cn Orbital Burns 

CD 

o 

RENDI 

VPG 

^ Time of 1st Burn 

CD 

O 

p— 

w Time of Last Burn 

CD 

O 

r—i 

d 

-P 

•rH 

rQ 

P 

O 

O rj 

> tf 
<3 CQ 

(ft/sec) 

Note (See Comments) 

a> Time of First Burn 
o 

w Time of Second Burn 

CD 

O 

cn Time of Last Burn 

CD 

o 

r*H 

d 

P 

• H 

rH 

0 

p 

o “ 

< -M tf 

<J pq 

(ft/ sec) 

Note (See Comments) 








4032 

25050 

17697 

(i) 

38553 

58382 

102716 

14138 


. i 

0 

90 

28. 56 

483 

100 

14,444 









HI 








43488 

75863 

15489 

(2) 

6846 

25810 

64800 

14032 

ESI 

2 

14502 

74. 92 

0 

15089 

100 

8762 





19103 

27383 

30125 

6057 

(5) 


(6) 






(7) 


(7) 







3 

19260 

84. 84 

3. 05 

19764 

276 

10127 

21372 

-- 

7378 

_1 


L 






RENDI 

(1) First, Min- AV Found 

(2) Best Min- AV Found 

VPG , 

(3) Bielliptic Lob 

(4) Full Orbit Phasing Chase 


(5) VPG Full Orbit Phasing Lob 

(6) Launch site is 80. 25 degrees from 
Target's Ascending Node at Launch 

(7) One OrbitaL Burn At Target Altitude 
to Match Velocities. 


> 

hH 
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Determine Target Ephemeris 
at Reference Time (to) 

Call EPHEM (R J0 V T0 ', t T0 ' )/x , 
DUM1, to, 0.001, 1 . , _r, _v, TOF, 
ma, a^} 



© 


Note:' Numbers on Blocks Correspond 
> to Statement Numbers in 
Fortran Source Program 
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A 


[ Calculate Times to Launch with Minimum 

out of Plane Angle (r 

MAX ' = MAX 

^ T 2 ) 

MIN ' - AL MIN 


Gall Launch (-1, A L , 

y, i T » t > to. 

MAX '• - A - L MIN '» T 1 

, t 2 j A l , Ai , A 0, 

• Ul > Ul) 





t’ = 




1 2 ' : 




\ 

r 

3022 

Do this block I = 

1, N 


Call Launch (1, \ L 

» ^LO » 

^T» W* *(!)> 


to, AL m ax 1 » MIN 1 j T t » T 2 > -A-L * A i , 

A9,_r L ,_v L , _n L ) 
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3028 


Following Selects I From Block 3022 which yields Min. out of 
Plane Error ( Ai ) and Secondly Minimum Launch Time (t) 


KNT =0 



3038 



Min= 1 '3~<IAi{l)|-lAi{2}j>—Min -2 


lAi(i)MAi(3)i: 


KNT = KNT + 1 j Min = 3 
Min Ait (KNT) =3 KNT =0 


~<lAi(l)|-|Ai(4)l>- 


3044 


Min =4 
KNT =0 


3052 


KNT = KNT + 1 
Min Alt (KNT) = 2 


3048 

\or 


N - 3 


3050 


KNT = KNT + 1 <3— / 

Min Alt (KNT) =4 #<!Ai(3)|~IAi(4)|>2 


Min =4 





r<TAi(2)I-lAi(3)i 


Min = 3 
KNT = 0 


N - 3 


KNT = KNT + i 
Min Alt (KNT) = 3 

3062 


N - 3 


3072 


lAi(3)UAi{4)|: 


3076 
Min = 4 


1Ai(2)MAi(4)|: 


_ 3074 | 

KNT = KNT + 1 
Min Alt (KNT) = 4 


3068 



Min = 4 
KNT =0 


l 3100 ) 
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^Selects, 

Minimum 

Launch 


Calculate Launch Conditions 
at Time t 

Call Launch (1, ^ LO > ^ t, 

V AL MAX’ AL MIN’ T l* T 2* AL ’ Ai> -L’ 
—L* ^L } - 


DO 3110 I a 1, KNT 
J S MIN ALT (J) 



MIN a MIN ALT (I) 


3110 


t*t (MIN)* 

A = A (MIN)* 

L L 

A i = Ai (MIN)* 
A0 =A0(MIN)* 
t = t. 


Write Out 

£j, A L> Ai, A@ 


R , = R r 
-l —E—I 


-1 " -L 

r a R 
a o 


r a R 

P o 


& T a fl + W (t - t ) 

L LO E ' 1 O' 


Calculate Ascent to Parking 
Orbit Parameters 


Call UP1 ($100) 
*2 = ‘l + AT 



5z = £ 


-Z = z 

A V = A V. 


Write Cut 

*2- Av „- A 'V a 2 . y 2 


Initialize For Rendezvous 
Intercept Calculation 

CaU EPHEM (R^, V^, t<> , M , 

A0, t 2> 0.001, r , v , TOF, ma. 


Call RENDI (R ? , V ? , r, v^ 0.., 

3 ' °- V A V U .H) 
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t = t l 


Call Launch (1, X , , Cl , i^t 

fo, al ma} , t , 

Zis -A/ ~L ] 

A I, t.j.T-, AL, 

MIiN , 1 2 

_r f a COS AS 

r + SIN A (9 n 

8 f ' a COS' 1 

[L •^TC ) /|R TO |]* 

SIGN £ 

®-TO 


Call EPHEM (R to , V^, 

fi,- 6 f ', DUM, 0.001, 0, R. Tf , 

_V , AT, ma, DUM) 

t* = t +AT 

f o 

-Sf = -St£ 

Af Tf 


CED 


AS, 


1-2 Main Program 
heet 7 of 1 2) 
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Intercept No 
*Mis sion 


V f = 0 


2032 


Last Av\ 
Remaining Less 
sThan Previous, 
Low 


Call UP2 ($100] 


mam 

■R 

= r 

-*■*2 


V- 

= v 

—2 


t 

= t, +At 

2 

1 

R f 

= r 




= v 



t 

= t/ 

f 

f 


Write Out 

m, tj, . f , a l , a V r , e c 


t* ~ t 

a l* = a L 

Ay* = Iav 


a l = a l 


AL ' 
MAX 


2038 
+ Aa l1 

2040 

1 = al| 


AL 1 =AL 
Mj.N 


0, + 

j 2048 
t l =t l +At l 


2050 


nxceea 

Allowable 



t ! ^ fc * 


Remaining 
Negative « 


Launch ^ 
Angle Exceed 
Maximum > 
^Allowable 


A = A 
L L MIN 


LC1 = -1 


Write Cut 
No Solution 
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^Surveillance] 



L MIN“ 

*( 7T- A )• 

- L MAX ' 


0, + 


4004 


a l* ~ A L min 

Lib 


1 -- £*2— 



L MAX 


4006 


i, • = COS ” 1 1 

MAX l 

CCSX I SIN A l *] 




4008 

■|Targ. Lat|v. No 



5000 ^ 


Calculate Parking Orbit 
M = 4 

A s A # 

^ T T 

a l nm - AL 
a L max’ = al 

l T " X MAX 
= 0 

Call Launch (1,XL, XiLO, 

jAL, Ai, A0, r , v , n ) 
jU 
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5000 



n 

= -n 



= 7r 

w 


£2 T = Sign (J^ • ji) tcos 


-I r K X n 
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t = t. 


Call Launch {1 # 0, XL, iQ> r 

, JL/W JL 

1 T ’ t0> AL MAX’ AL MIN’ T l s r 2* AL ’ 


Ai, A 6, _r L , _v , ja, ) 

r = R 
a o 

r - Rrp 

P T 

n dn = ( -C°Sa r! -SINi2 T , 0) 

‘ 6 = COS' 1 (r SIGN fin 

p —l Dim l^- 

■Call UP1 ($110) 

~1 5038 


•DM X -T- 


= r 

= V 


= v 

t, = t +At 

Av =Av 

O 

Av =Av d 

n R 


5056 


r 

A * = AL 

No 

L 

3- 

8* = 8 

cn 

O 

O 

00 

i . .= i 


MAX T 


cos A n= 


^8^/Yes _ 

5064 


(cos 8 -sin 2 xJ 


cos 2 X. 


An= cos ' 1 [cosAXl ] 

At £ = (4 a -) 4SIaN (t 3 - V 

» TT ' 


v E 7 
tf = t f +At f 


"AvV Yes 


»>-# t f = t f + 600 


5068 


-1 /—2 * -tN * r -s 

S V IVI. j SIGN [ £ V^r , --a] 

Call EPHEM (R ? , _V 2 , t , DUM1, 0.001, 

0,_r, _v. At, ma, DUM) 

t 3 =t 2 +AT - 

v , 3 *_v 

XI = £2 4-AU .+ w (t -1 l 
V.T3 LO T E l 3 , o' 

r T3 = cor Xj, CCpQ T3 , smfl , einX. 

8= cos-1 [(Rj •i r3 )/|R 3 |] 


Write 

No First Pass 
Solution 



Write Out 

V V AV o 

V a l 

*2» -&2* %Z‘ Av , 

V k 3 , X 3 , 8 


lSj 












AI-55 



Figure 1-3 Launch Subroutine 
(Sheet 1 of 3) 

















Figure 1-3 Launch Subroutine 
(Sheet 3 of 3} 






Figure 1-4 EPHEM Subroutine 
(Sheet 1 of 3) 











Figure 1-4 EPI-IEM Subroutine 
(Sheet 2 of 3) 











Figure 1-4 EPHEM Subroutir, 
(Sheet 3 of 3) 








Figure 1-5 RVVR Subroutine 
(Sheet 1 of 3) 











I 



Figure 1-5 RVVR Subrout 
(Sheet 2 of 3) 












AI-63 


£ 



r 



A-fL. - cos 6 


Test orbit to determine if 

j. a 

5 o 


hyperbolic, parabolic, rectilinear. 

Bound = ; —-x -- 

sin 6 


reentry, or elliptic. 

a 




This section of flow 
changes the value of 
rif necessary to 

make it compatible 
with the transfer. 


No 


/- 

l Return J 



.- . ▼ - 

X MASAV 


No 

.. 

L ~ i 



^ MASAY 

== B ound 4 1 



= Bound - I 


- cos e 

r a 


m = 


a sin 9 


+ X 


MASAY 


X = m 

MAS A / a 


Write "Slope Change to. 
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Input: 0, m ,r ,C,m 

O O cl 

Outfrat: tof 


Initialize 
NS = 5 ! 
tof = 0 


FOR ri = 0, NS 

cos n A 6 = cos (n-l}A0co.sA#-sin(n-l)A0sin A 9 
sin n A £? = sin(n-l)A#cos A 6 + cos (n-l)A0sinA@ 




) cos nA 9 — - ■ sin nA0 j 


AS r, 2 V NS * + 14 T‘‘ * T , * 4 2 , 

TT .l 14 L "n +16 £ , , r n - 7(r o Tr 2Ns' 

n = 0, 2, 4 n = 1,3, 5 


tof = El/C 
N 


+ 2 A0 (m r" 2 - m m ) j 
o a a 2NS J 


Itof - tof 


tof = tof 

N 

NS = 2 NS 

KOWNT = 
KOWNT + 1 


^ Is ^ 

KOWNT<3 


Return 


Write "Simpson’s 
Time Calculation 
off bv " 


Return 
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Inputs: _r i ,_r* i »_r T .»r f ^i» TSTOP, 
m .kj.NOIT.TSTART, AV M , p. 


Al-66 


Read in 

T , TRIM 
A 


Initialize: N180 = 0, IPL = 1 
TMINSA = 0,DLE = 10" 3 , XK=10- 6 , 
NN = 0, TO - 0, t = 0,p = 0 


Compute interceptor and target 
periods* Initialize AT by 
setting it equal to 1/4 the 
smaller of the two periods 


, 1000 


Call EPHEM subroutine Option B 
to compute the interceptor ephemeris 
V r 0 an< ^ target ephemeris 

for t sec after epoch (t = 0) 


Call EPHEM to compute target 

position and velocity r , r,* , 

r x a. 

and slope m A ^ after transfer 
A0 

taking T A sec, 

° A 


\ 

f 

Compute Ai: 

^a' =r TA + e 

ti Pointy 
r TA 7 i r TAl 

. 1 _ 


Call RVYR with slope m. 

A0 

to compute time ef flight T 


of interceptor from r to r 

—a 


m Al =m A0 +0 * 01 
KOUNT = 0 


Is Yes 

k ,< 0 


Compute interceptor position 
projection on target plane, r. 


Compute angle 9 between 

r. and r_ (0<6<2-rr) 

— inp —T v “ 1 



r 971 


Call RVYR with slope m , 

1 A1 

to compute time of flight T 

A. 

of interceptor from r to r 

—» o —► a 


Use regula falsi iterations to 

find zero of VALUE = T - T , 

A AI 

with m ^ as independent 

variable. Compute new value of 

xn. called m . . 

TA A2 


KOUNT - KOUNT + 1 

















m 


DLIM = 0. 5 - 8 


9< 2tt+ 0. 5 


Yes 

__1_ 

. . . p 

DLIM = 0 


t t+ lz 


JN 

1 

\ 

r<C 

L_ 

t + TSTART 
S^TSTOP/ 

Return 

1 |Yes 


ULIM = 2tt -0. 5 -0 


KOUNT>20 


Write : 

"Iterations did 
not Converge" 


Return 

No 


S' Is 
I VALUE | 
< T LIM 


6 ~ DLIM 
o 

= ULIM 


Call EPHEM Option B to com¬ 
pute the velocity r' of inter- 

—»• o 

ceptor at arrival point r 


D LIM < ULIM 


Call EPHEM Option A to .compute 
target position and velocity ,’ 5 

_iTA*-^TA slope m T A.* and time 

of flight T_ A _ resulting from 
° TAG 

angular transfer of DLIM from r,l 


r = r„, . 
a TA 


+p ' TA /|' 


Compute AY required at first 
and last burns and also total AV 


Call BEST to enter information 
required for finding minimum 

AV and compute new AT 


TA * TA 1 


<b 
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Call R WR to compute interceptor 

time of flight T _ from r to 
° AO o 

the aim point r 


Call RVVR to compute interceptor 
time of flight T and required 

velocity V from r to aim 
—r o 

point r 


KOUNT s= 0 
VALO ~ T. 


VALUE = T TA , -T 


/is\ 

Value 

voy 


VAUP = VALO 
THUP = DLIM 


VADOWN = VALO 
THDOWN = ULIM 


KOUNT = KOUNT - 1 



IOUNT 

o30 


VAUP=VALUE 

THUP-'TKl 


VADOWN = VALUE 

TKDOWN = TH1 

----- 



Test to see if 
23 zero occurs 


Write "Iterations 
did not converge" 

Return 


VALO .VALUE 

N \>o 

_ (No 

20 X 
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Input: DY, T, TH, IPL 
DT, T MINS A, NOPT 
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DVSAV = DVEL{I) 
DVEL(I) =DVEL(I+ 1) 
DVEL(I + I) - DVSAV 
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DVEL(I) = DVEL(I+ 1) 
TM(I) = TM(I + 11 


TM(I) - TM1 


1ST 1 


DENOM 


~ 2 ^DVEL(l) (t 2 -t 3 ) + DVEL(2) + DVEL(3) ( tj - r ? )| 


DENOM 


(r - r W r - r } ( t. - r, 
1 1 2V 2 3> v 3 1 


DVEL(l) (r/-^ 2 ) + DVEL(2) + DVEL(3) (r^-r^) 

DENOM 


T = TM(1) +T-ISTI T=TM{3 )+I3 TlSGN TM(3)-TM(2 j]| 


IPL=4 


I T~TM(3)| 


Return 


ino T = TM(3) + DT*SGN|T-TM(3)J 


IPL = 4 


Return 
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Figure I-10 VPG Subroutine 
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Figure X-il PCA Subroutine - H 
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Write: "Does 
not Converge” 
NMM = 1 


Return 


Fig 
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Appendix II 

TOP - TRAJECTORY OPTIMIZER PROGRAM 

2. 1 INTRODUCTION 

The Trajectory Optimizer Program (TOP) is a digital computer 
program written in FORTRAN IV to run on IBM 7090/7094 computers. The 
program is used in mission analysis and planning to generate local minimum 
time or minimum fuel trajectories from launch through various orbital ma¬ 
neuvers. 

The program employs a' direct search optimization algorithm 
which systematically varies the paramater of an initial (input) trajectory in 
order to minimize a pay-off function - the time or AY required to complete the 
mission. Inequality constraints on trajectory parameters are enforced in the 
search routine while functional constraints are handled by the addition of 
penalty terms in the calculation of the pay-off function. 

This program was written to demonstrate feasibility of the Tra¬ 
jectory Optimization Mode of the On-Board Mission Planning Function (see 
Section V of the main report). 

2. 2 PURPOSE AND USE 

TOP represents a flexible and efficient technique for generating 
local optimum time or fuel mission trajectories. The trajectories to be 
optimized can start at launch or from orbit and can include any number of 
orbital burns for orbit transfer with or without rendezvous (the present version 
of TOP is limited to four orbital burns). All orbital burns are treated as 
impulses while the ascent-to-orbit employs an explicit guidance algorithm 
which provides nearly optimal performance. 

This program is used to optimize the trajectories generated by 
the Trajectory Generator Program (Appendix I), but it can be used with any 
initial trajectory if the proper trajectory parameters are) specified (see Para¬ 
graph 2.5 for input requirements). ~ ; r 
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2. 3 PROGRAM FEATURES 

The present version of TOP does not represent a production level 

9 

program. It does represent a possible technique for on-board generation of 
locally optimal mission trajectories. The program is the result of different 
efforts of Task II during this phase. The present program configuration is 
not designed for high efficiency, and certain approximation used in it must 
be studied in more detail before a final version is constructed. 

Some of the features of the present version are: 

© Ascent to orbit or intercept (e.g., direct-ascent 
rendezvous capability). 

• Spherical gravity model. 

• Kepler arcs for coast phases. 

m Total AY or total mission time can be minimized. 

• Maximum AV constraints enforced by a penalty term in 
the pay-off function (OBJECT) when minimizing total 
time. 

• Components of the argument vector (U) can be con¬ 
strained in the search routine (TEXT). 

• Up to four orbital burns can be considered, 

• Transfers can be with or without rendezvous. 

m Near-optimal vacuum ascent trajectories generated by 
an explicit guidance algorithm. 

• Atmospheric ascent trajectories calculated with 

analytical equations and vehicle-dependent empirical 
parameters. 

• All orbital burns treated as impulses. 

• Initial parking orbit, after the ascent phase, is circular. 

• Impulse splitting cases can be handled. 

• 180 degree transfers are allowed. 

• Optimum transfer plane determined for 180 degree 

transfers. 
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Some of the approximations implicit in these features are discussed in Paragraph 
2.4. 3 of this Appendix. 

9 

2.4 PROGRAM OPERATION 

The over-all operation of TOP is illustrated in Figure II-1. The 
program consists of a Main Program (MAIN) and two major Subroutines, 

TEXT and OBJECT. The Main Program reads in all data, performs neces¬ 
sary conversions, sets up program paths, calls OBJECT and TEXT, and 
formats all printout. TEXT is the optimization subroutine which itself con¬ 
sists of a number of subroutines that cam be used in various combinations. 
OBJECT is the subroutine which calculates the value of the pay-off function 
(P) corresponding to the present value of the trajectory parameters (CJ) as 
determined in TEXT. 

Referring to Figure II-1, MAIN reads in all data including the 

initial trajectory, which is specified by a set of parameters (TJ (I) , I = 1, . . . 15). 

The initial trajectory specification, U^, is used with OBJECT to provide the 

initial value of the pay-off function, P . TEXT is then called to begin the 

o 

search procedure, and control remains in TEXT until a solution has been 
found. At this time, MAIN formats and prints the output data and either 
proceeds to the next case or stops. In TEXT, the OBJECT subroutine is 
called to provide a new value of P whenever U has been changed. 

2.4.1 MAIN 

A math flow of the executive program MAIN-is presented in 
Figure II-2. The symbols employed may be found in Table II-1. In 
addition to reading input, printing output, and calling OBJECT and TEXT, 
this program: 

• Sets constraints on appropriate terms of the argument 

vector (U) with the read-in values of U (I) and U . (I) . 

— max mm 

• Identifies the components or U to be searched over 

: - '* (BVA(I) = 1 if U (I) is to be varied, BYA (J) = 0 if U (J) 


is to be held fixed). 
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• If an impulse-splitting maneuver is being considered, 

MAIN reverses the order of the orbital burns before 
calling TEXT and changes them back after returning 
from TEXT. 

The last item relates to the problem of determining the velocity required to 
transfer between two points in space in a finite time when the points are 
coincident. This is the case for the so-called impulse splitting maneuver 
which is employed in rendezvous missions to introduce a phasing orbit. 

When OBJECT is first entered to determine the initial value of the pay-off 
function a check is made to see if the last two burns are within p (RHO)-feet 
(an input) of each other. If they are. an indicator is set which, on the return 
to MAIN,. causes the order of the orbital bui&ns to be reversed-(i-. e. , the 
program starts at the target orbit and works back to the parking orbit). This 
is shown on Figure II-2 just before TEXT is called. The optimized orbital 
burns are then changed back to the proper order before final printout. By 
reversing the order, the first two orbital burns are now coincident in space 
but, because of the technique employed in OBJECT to propagate trajectories, 
this causes no problems. 

2.4.2 TEXT 

TEXT is the optimization algorithm subroutine. It consists of 
several options which can be used in various combinations. This subroutine 
is completely described in Appendix IVwhich includes Math flows, Input/ 
Output, and examples. In Appendix IVthis algorithm is entitled DSOP 
(Direct Search Optimization Program) and its various options are described. 
Only the Pattern Move Search (PMS) and IJNIVAR options have been em¬ 
ployed in TOP because of the ability of this combination to handle argument 
constraints. Any search type algorithm could be employed with TOP but some 
reprogramming would be required to interface it with MAIN and OBJECT. 

2. 4. 3 OBJECT 

The payoff function associatedly/ith a/parthcnilar trajectory is 
computed in the OBJECT subroutine. The input to this routine is the 
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argument vector U which describes the trajectory and the output is the corre¬ 
sponding value of the pay-off function P. The argument vector is supplied by 
TEXT while OBJECT calculates the time or AY (the presently available pay¬ 
off options) required on that trajectory. 

OBJECT makes use of five subroutines: LAUNCH, UP1, UP2, 
EPHEM, and RVVR. The launch conditions are calculated in LAUNCH as a 
function of the launch site and launch time. UP1 and UP2, employing the 
explicit guidance algorithm described in Par. 6.4 of the main report, 
generate near optimal ascent trajectories to orbit injection or intercept. 

The determination of position, velocity and time on coast arcs performed in 
EPHEM. The RVVR routine is used to calculate the velocity required for 
the orbital burns. These subroutines are also used in the Trajectory 
Generator Program (TGP) and a complete description of them and TGP may 
be found in Appendix I. 

Figure II-3 is a math flow of OBJECT and the symbols are 
defined in Table II-1. 

The OBJECT subroutine employs certain approximations in order 
to provide an efficient technique for generating ascent and orbital segments 
(this routine is called hundreds of times in the process of optimizing a 
mission trajectory). These major approximations are related to three areas 
of trajectory generation: Ascent, Coast, and Orbital Burns. 

• Ascent - the ascent segments are generated with UP1 or 
UP2 depending upon whether orbit injection or intercept is 
required. Both of these routines (which are really two 
entry points to one program) use the explicit guidance 
equations, developed during this study, for the vacuum 
phase and analytical equations for the atmospheric phase. 
The form of the explicit equations are such that they 
introduce almost no performance degradation when com¬ 
pared to an optimal solution. The analytical equations, 

: w which are necessary to initialize the explicit guidance 

equations, are merely spherical trigonometry equations 
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relating conditions at first (zero) stage burnout to the 
launch conditions. The equations require three vehicle- 
dependent parameters which are used for all launch con¬ 
ditions. This analytical technique of handling the atmos- 
. pheric phase of ascent is very efficient computationally but 

it does introduce prediction errors of the order of 1 - 3 
seconds in the injection time. This aspect of the ascent 
phase can and should be improved with additional study of 
efficient techniques for generating conditions at the !! top 
of the atmsophere ,f as a function of launch conditions and 
vehicle configuration. 

• Coast - All coast arcs are computed with an ephemeris 
routine (EPHEM) which employs Kepler-type equations, 
that is, assumes a spherical earth. In the case of long 
mission times, this will result in sizable errors if un¬ 
corrected. The present program could be upgraded by 
considering the effects of oblateness on coast trajectories. 
Closed-form expressions for these effects are available. 

• Orbital Burns - All orbital burns after ascent are treated 
as impulses. This approximation is justified because of 
the almost negligible performance penalty incurred. The 
results of Reference (1) , for example, indicate this 
penalty is less than 1% for an extreme case (AV = 10, 000 
fps. ) and provide a prescription for generating a nearly 
optimal finite-thrust trajectory for data derived from the 
impulsive trajectory. 

One other feature of OBJECT should be mentioned. This involves 
the program's ability to generate 180-degree transfers when the transfer plane 
is undefined. This would be the case, for example, for a simple Hohnann 
transfer between non-coplanar circular orbits. Specifying the times 
(positions) of the two burns is not sufficient because the transfer: plane is l 
still arbitrary (i. e. , the plane change can all be made at the first burn, or 
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all at the second burn, or split between the two burns in an infinite number of 
ways). The logical solution, of course, is to^divide the plane-change between 
the two burns so as to minimize the total AV required for the two burns. This 
is exactly what is done in OBJECT when the trajectory specification, U vector, 
calls for two burns that are within XK1 -degrees (a program constant, presently 
set at approximately 0. 08°) of being 180 degrees apart. When this situation 
arises, a simple one-dimensional search is performed over the angle de¬ 
scribing the transfer plane, to find the value which minimizes the AV required 
for these two burns. This capability is required in order to accept such trans¬ 
fers as part of an initial trajectory specification. The VPG Bielliptic option, 
for example, always includes such a transfer. 

2. 5 PROBLEM SPECIFICATION 

The coordinate system used in defining input and output mission 

data, and in which the problem is solved, is the XYZ ECI (Earth Centered 

Inertial) system illustrated in Figure II-4. The XY plane is the equatorial 

plane and the Z axis is through the north pole. 

If the initial trajectory starts at launch (LCL = +1) the launch site 

is specified by four parameters which are varied during the optimization. 

t = Time of launch = U(12) 

L 

A = Launch Azimuth = U(13) 

L 

R = C ircular parking orbit radius = U(14) 

P 

a = Dog-leg angle - U(15) 

L 

The present program only handles circular parking orbits. The dog-leg 

angle a denotes a rotation about the launch radius which defines the de- 

sired ascent plane after the atmospheric phase. A defines the plane for the 

L 

atmospheric phase. 

The target orbit is specified by the target 1 s ephemeris at some 

epoch: 

= Target 1 s position at t T 
o o 

V = Targets velocity at t T r; 

T o u 
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trp - Target* s epoch 

o 

If the initial trajectory starts from a parking orbit, LCL = -1, this orbit is 
defined by appropriate ephemeris data: 

= Spacecraft position at t^ 

= Spacecraft velocity at t^ 
t = Spacecraft* s epoch 

The specification of the orbital burns requires 3+4 (N-2) para¬ 
meters, where N is number of burns (N >2). If rendezvous is required, one 
less parameter is involved because of the phasing constraint (i.e. , the space¬ 
craft and target must be coincident at the time of the last burn). For all 
values of N, the first and last burns are defined by: 

t^ = Time of first orbital burn = U(l) 
t^ = Time of last orbital burn = U(2) 
ma^ = Slope ox arrival at last burn = U{3) 

The position of the first burn is determined by t^ and the initial parking orbit 
(input or result of the ascent phase). The position of the last burn is de¬ 
termined by t^ and the target*s ephemeris. For a two burn orbit transfer, 
this and the slope (tangent of the flight path angle) at arrival at the last burn 
are sufficient to determine the two impulsive burns required. In the case of 
rendezvous, ma^ is not used; the phasing constraint is employed instead. 

When there are more than two orbital burns, four more para¬ 
meters are required for each additional burn. The parameters used are: 

AV^ = Velocity impulse at i^ burn = U(4i), U(4i+1), U(4i+2) 

A£h = Range angle between i** 1 and (i+l) st burns = U(4i+3) 
where i = 1, 2, ••• (N-2). This combination allows for efficient generation of 
the intermediate trajectories, because only the ephemeris routine (EPHEM) 
is required. ' For example, t and the spacecraft’s ephemeris with EPHEM 
yield the position and velocity just before the first impulse (R^ and V ). 

Adding AV^ to gives the velocity after the impulse which, with R- and 
A0^ , is sufficient for EPHEM to generate the position and velocity just 
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before the second burn (R and V ). Adding AY to V , proceeding as before 
enables the next arc to be computed, and so on. 

The initial trajectory (input) is thus specified by assigning the 
proper values to the initial value of the argument vector U(i), i = 1, 2, « • • 15. 
In addition, other input data is required to completely specify the problem. 
The total input required is indicated in Table II-l. 

2.6 TEST CASES 

To demonstrate the capability and performance of TOP, four test 
cases are included. Table II-2 contains a description of each case, and 
Table II-3 illustrates the results obtained from TOP. 

Case 1 - This is a two burn orbit transfer v/ithout rendezvous 
between two similar orbits inclined at 5 degrees. The initial (input) tra¬ 
jectory corresponds to a perigee-to-perigee transfer arriving at the final 
orbit with zero flight path angle. This maneuver requires 7619 ft/sec of AV. 
The optimized trajectory requires 5663 ft/sec, a saving of 1856 ft/sec. This 
case required 155 evaluations of the pay-off function (i.e. , passes through 
OBJECT) and 2. 78 seconds of 7094 time. The running time figure contains 
I/O operations for writing input data, intermediate results, and final results. 
This case corresponds to.Optimum 4, Table 1, page 1868 of Reference (2) 
which lists a value of 5800 for AV. This difference in results is due to the 
fact that the values for AV in Reference (2) are obtained by reading 
contour plots which have a contour interval of 500 ft/sec. 

Case 2 - The problem here is to rendezvous with a target in a 
highly elliptical orbit with a large semi-major axis. The initial orbit is a 
100 n, m. circular parking orbit inclined at 30 degrees to the target orbit. 

The starting solution is a three-burn bielliptic transfer which was generated 
with the Variable Point Guidance (VPG) routine. This routine is used in the 
Trajectory Generator Program and is described in Appendix L The VPG 
planned trajectory requires 13, 528 ft/sec of AV as compared to the optimized 
results of 12,715 ft/sec. The very large number of function evaluations., 
(1090) are a result of slow convergence in the algorithm and, especially, 
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because during many of the evaluations 180 degree transfers are involved and 
the optimum transfer plane must be determined (see Par. 2.4. 3). 

Case 3 - This is exactly the same as Case 2 except that the 
starting solution is a three-burn, full orbit phasing result from the VPG 
routine. In this case the input contains an impulse splitting maneuver be¬ 
cause the second and thrid burns occur at the same point in space. The 
technique described in 2.4. 1 was employed to handle this case. The opti¬ 
mized trajectory requires 11,496 ft/sec as compared to 11, 764 for the initial 
trajectory. This case, in particular, demonstrates the near-optimal per¬ 
formance of VPG in certain applications 

Case 4 - The target orbit in this case is the same as Cases 2 and 
3 but now the missions start at launch. The input trajectory was generated 
with the Trajectory Generator Program (Appendix I) and involves a launch 
at t = 0, at an azimuth of 90 degrees into a 100 n. m. circular parking orbit 
with no dog-legging (U(15) = 0). This ascent trajectory results in the minimum 
out-of-plane angle with the target orbit that can be obtained from the given 
launch site (latitude of 30 degrees). The input specification then involves a 
two-burn rendezvous maneuver with the target which results in a total AV 
of 31, 673 ft/sec and requires 15, 706 seconds of time. The optimized tra¬ 
jectory required 255 evaluations of the pay-off function (total AV) and re¬ 
sulted in a mission requiring 14, 769 seconds and 31,417 ft/sec of AV. This 
is a 256 ft/sec savings in AV and 937 seconds in time. The optimized mission 
involves a launch at t = 0 (launch time was constrained to be positive) at an 
azimuth of 92. 7 degrees with a 1.287 degree dog-leg after the atmospheric 
phase. The parking orbit altitude was reduced to its constrained minimum 
value of approximately 90 n. m. The two orbital burns were altered in time 
of occurrence, and thus in position because the parking orbit is different. 
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TABLE II-1 

TRAJECTORY OPTIMIZATION PROGRAM 
SYMBOLS 


Math * 
Symbol 

Description 

FORTRAN 

Program 

Symbol 

LCL 

I. Inputs to Main Program 

Equals +1 for Ascent, Equals -1 for maneuvers 

that start from Orbit 

Note: The following inputs are required only for 

LCL = +1, i. e. for ascent. 

LC L 


Launch site latitude (degrees) 

LAML 

r Lo 

Launch site longitude relative to x-axis at t 

LU 

(degrees) 

OLO 

t LO 

Reference time 

TLO 


Longitude of ascending node of target orbit (degrees) 

OT 

AD, _ 

T 

Longitude of target overflight point relative to 

launch site (degrees) 

DELOT 

X -p 

Latitude of target overflight point (deg. ) 

LAMT 

*T 

Inclination of target orbit plane to x - y plane (i> 0) 

EYET 

*> 

Radius of target overflight point (ft. ) 

RTP 

AV 

max 

Maximum AV available (ft/sec) 

DYM 

‘l 

Time of launch 

U(12) 

a l 

Launch azimuth (degrees) 

U(13) 


c 
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Table II—1 Trajectory Optimization Program Symbols (Continued) 


Math * 
Symbol 

9 

Description 

FORTRAN 

Program 

Symbol 

r p 

Parking orbit radius {ft. ) 

U (14) 

“L 

Dog-leg angle (degrees) 

U(15) 

M m 
. T 

Equals 1 for Intercept 

2 for Direct Ascent Rendezvous 

3 for Orbital Rendezvous 

4 for Reconnaissance with parking orbit 

5 for Reconnaissance with first pass 

overflight 

6 for Orbit injection using a parking orbit 

MT 

hmin 

Lower bound on t 

L 

AC(i) 

^Lmin 

Lower bound on A (deg. ) 

L 

AC(2) 

^Lmax 

Upper bound on A (deg. ) 

L 

AC(3) 

R Pmin 

Lower bound on Rp (ft. ) 

AC(4) 

R Pmax 

Upper bound on Rp (ft. ) 

AC(5) 


Switch to control the use of above bounds. 

Note: The remaining inputs pertain to orbit maneuvers. 

NAC(1 - 5) 

N BURNS 

Number of orbital burns (0, 1, 2, 3, or 4) 

NBURNS 


Desired orbital maneuver. Equals +1 for 

rendezvous and -1 for orbit transfer. 

MANEUV 


Desired minimum Equals +1 for minimum fuel 

and -1 for minimum time. 

MINREQ 

r 

o 

Spacecraft position at time t > 

RO 
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Table II-1 Trajectory Optimization Program Symbols (Continued) 


Math * 
Symbol 


Description 


FORTRAN 

Program 

Symbol 


Spacecraft velocity at time 


Reference time for r and v 

— o —o 


Target position at time t r 


RTO 


jTo 


Target velocity at time t r 


YTO 


Reference time for r and v m 

-»T.o _ 4 ,To 


TTO 


Time of first orbital burn 


Time of last orbital burn 


U(2) 


Slope at last orbital burn 


U (3) 


Imoulse at first orbital burn 


U(4), U(5), U(6) 


Range angle between first and second orbital burns 


AY 

-» 2 


Impulse at second orbital burn 


U(8), U(9), U(10) 


Range angle between second and third orbital burns 


U(ll) 


f max 


Upper bound on t 


AC(6) 


Switch to control use of AC(6) 


NAC(6) 


AV max 


Maximum available AV 


FC(1) 


Weighting factor to control use of FC(1) 


NFC(l) 


Gravitation constant 


U, GC 


Rendezvous iteration time tolerance (1 sec. ) 


EPT 


Impulse splitting indicator (5000 ft. ) 


RHO 


Optimization search tolerance ( = 1) 


BIAS 


End of inputs 
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Table II—1 Trajectory Optimization Program Symbols (Continued) 


Math * 
Symbol 


Description 


FORTRAN 

Program 

Symbol 


II. Internal and Output Symbols to Control and 
Object Programs 

-5 

Rotation rate of earth = 7. 292115- 10 radians/sec WE 


NPASS 


.i=l, 4 


. i = 1, 4 

l 


T i = 1, 4 


*5' J 5 ' T 5 


Earth equatorial radius 

Range angle of S/C from epoch to first burn 
position 

Range angle of target from epoch to final burn 
position 

Total velocity change 

Number of presses through Object Subroutine 

th 

Position vector of i orbital burn (ft. ) 

th 

Velocity vector of i orbital burn (ft/sec) 


Time of i orbital burn (sec) 


Position, velocity, and time at cut-off for 
ascent 


THF 


NPASS 


R( , i) 


V( , i) 


T (i) 


:( ,5), V( ,5) 


|AVp| i = 1, 4 Magnitude of velocity change due to i orbital 

burn (ft/sec) 


DELV ( , i) 


IAV I 

_•» D 


Magnitude of velocity change needed during ascent DELV ( , 5) 

(ft/sec) 

Out of plane angle between S/C orbit (or launch plane) DELI 

and target orbit (radians) 
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Table II— 1 Trajectory Optimization Program Symbols (Continued) 


Math * 
Symbol 

Description 

FORTRAN 

Program 

Symbol 

a e 

Range angle from launch site to intersection 

with target orbit (radians) 

DELTH 

r T 

Unit vector through launch site 

RL 

—*-L 


Jl 

Unit vector in launch plane and perpendicular 
to r^ 

VL 

Jl 

Unit vector normal to launch plane such that 

n = r x v 
— » _L/ — ** Lj 

NL 

Jl 

Unit vector normal to orbital plane 

UNIT 1, N] 


6 P 

Longitude of perigee from descending equatorial 

node (radian^) 

THETP 

r T ^ 

Position vector of launch site 

RLB 

— LB 


r 

a 

Orbit apogee radius 

RA 

r 

Cut-off position vector (ft) 

RC 

c 



V 

Cut-off velocity vector (ft/sec) 

VC 

_«. c 



av r 

AV remaining at cut-off (ft/sec) 

DELVR 

AT 

Cut-off time relative to launch 

DELT 

TF 

Desired time of intercept relative to launch (sec) 

TF 

Jr 

Final intercept position (ft) 

RR 

Jv 

i 

S/C velocity at intercept (ft/sec) 

VV 

JVT 

Target velocity at intercept (ft/sec) 

! 

i .. ! 

VVT 
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Table II-1 Trajectory Optimization Program Symbols (Continued) 


Math * 
Symbol 

Description 

FORTRAN 

Program 

Symbol 

a 

Semi-major axis of orbit (ft) 

A, AX 

TOF 

Time of flight (sec) 

TOF 


Range angle (radians) 

TUT 

VFIN 

Arrival velocity of S/C (ft/sec) 

VFIN 

m 

a 

Arrival slope 

MA 


( 


* In some cases tilex'*e is no math symbol corresponding 
to a FORTRAN program symbol. Such cases arise 
when in the math flow chart, descriptive phases are 
used instead of FORTRAN symbols. 
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TABLE II-2 


TEST CASES 


Case 

Initial 

State 

Final Orbit 

Out-of- Plane 
Angle (deg. ) 

Comments 

1 

Elliptical Parking Orbit 

Elliptical Orbit 

5. 0 

Two-Burn Orbit 


Semi-Major Axis = 5208 miles 

Semi-Major Axis = 6250 miles 


Transfer Without 


Eccentricity = 0.2 

Eccentricity ~ 0.2 


R endezvous 

2 

100 n.m. circular 

parking orbit 

Semi-Major Axis = 8444 n.m. 

Eccentricity = 0.5 

Line of Apsides 30° from 

Line of Nodes 

30. 0 

Three Burn Rendezvous 

with Target in Final 

Orbit 

3 

Same as 2 

Same as 2 

30. 0 

Same as above 

4 

Launch Site Latitude =30° 

! * ; i 

Same Orbit as 2 and 3, 

Lying in the Equatorial 

Plane 

30. 0 

(This is the 

minimum 

value obtain¬ 
able from the 

Launch Site - 

if Launch 

Azimuth = 90°) 

i 

Ascent to Parking Orbit 

and Rendezvous with 

Target in Final Orbit 
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TABLE II-3 


RESULTS OF TEST CASES 


1 Input 
Output 

2 Input 
Output 

2 Input 
Output 

4 Input 
Output 


rt * 

> £ 

.r-4 M 


Time c 

Burn - 

Time c 

Burn - 

Slope c 
Last B 

X - Co 
Orbital 

Y - Co 

0 rbital 

Z - Co 

Orbital 

U(1) 

(sec) 

U(2) 

(sec) 

U(3) 

U(4) U(5) U( 6 ) 

(ft/sec)(ft/sec) (ft/sec) 



U( 7 ) U(12) U(13) U(14) U(15)_J AV Passes Time 
(deg) (sec) (deg) (ft) (deg)[(ft/sec)_ 


2213 0 

2116 .07216 


7515 

19920 

- 

-173.9 

-4317 

-2108 

208. 7 

- 

- 

7825 

19114 

- 

2262 . 0 

-4797 

-645 

189. 9 


- 

2305 

33026 

- 

358. 4 

-5668 

-2788 

203. 2 

- 

- 

2370 

33373 

- 

879. 3 

-6329 

-1070 

201 .7 

- 

- 

3996.2 

15706 

- 

- 

- 

- 

- 

0 

90 . 0 

3970.2 

14769 

- 

- 

- 

~ 

- 

0 

92. 7 

Notes: 


7619 

5663 155 2.78 

13528 

12715 1090 

11764 

11496 234 9.83 

31673 


(1.) Not required for Rendezvous Cases E 
(2.) Used for Three-Burn Cases ^ 

(3.) Includes I/O Time to Write Output 


















Set for I-lj 15 

'2 O 

UMAX(I)=10 

38 

UMIN(I)=-10 . 


Doe s\ 
t£ have 
'Upper Bound 


UMAX(2) 

l fmax 


Ascent 


Dower Bound 



© 


UMIN(IZ) 

^Lmin 


UMIN{ 13) 

Dm in 
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UB(I)=U(I) 
1=1, 15* 
ISP=1 


Maneuver \.Q r1:> ^ Transfer 
Type ? 

TRendezvous 


Call 

OBJECT 


Ascent 


XI- XI + 

L LO 

"E X (t L‘ To* 


Convert fi 

j_j 

and X^ from 
radians to 


degrees 


< Im pul s e 

Splitting , 
Required? 


U(2)=U(1) 

U(1)=UB(2) 

U(4H 

u(5) > = av. 

U(6)J 

U(8) 

U(9) ^=AV 

U(10)J 

U(7)=2tt=A0 


NBURNS 


Resets 

Search 


NBURNS -1 Variables 
a for Text 


' ' 1-2 
U( 11)=UB( 11) 
UMAX (J )=UMAX (2) 
' 3 8 

UMAX(2)=10 

UB(I)=U(I) 

1=1, 15 


Impulse 

Splitting 

Required 


i=l, NBURNS 
R’.=R. 

— i -*i 

V!=V 

AV ! . = AV. 

T*.=T. 
i i 

| AVl^jAVI. 

R i _R 'NBURNS~i+l 
V V ! 

—i — NBURNS =i + l 
AV = AV 

— i — NBURNS-i + 1 
T i =T ’NB URNS-i+1 
l^liii^ll'NBURNS-i- 


Call Optimization 
Search Routine, TEXT 


Convert Sj, 8 f , AS ^ AS^, 
a fron radians to degrees- 


Output: t x , t f , m af , AV^, A0 

~Z* Ad Z-3* t L* A U’ R P’ “ L ;, ' R i* 
T.J AV !(i=l, NBURNS); 


•AV, 0 , 0 f ,NPASS 


Loop back and read new 
set of input cards 
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Maneuver 

^ T yP e ^ 


Rendezvous 

Intercept 


lOrbit Transfer 


Call RVVR 

o. X K1 »^» norb * 

TOF, v r , KKj 

, ) 
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(STOP) , 

Figure II-3 OBJECT Program 
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— NBURNS 

V -VFIN 

--NBURNS f. 


Compute 
|AV _ 


NBURNS 


Write 
"Iterations 
did not 
Converge" 


av k := 10 ' 


^ K A r --I K 

Compute | AV t 


— K 


NOP=0 


Call EPHEM Ootion A 

Sk'-Yk’-L* W’ 0 a’ 

T , AE, NOPj VFIN, 

O <r ^“ ; —«* 


A = |AV 


-M AV 


NBURNS 


^ transfer ^ 
been optimized 
with respect 
to a? 


Av = 

— NBURNS 

, V -VFIN 

— NBURNS U 


180° Transfer Case 
Optimization of AV 
with respect to 
plane change, a . 
Three successful 
passes used to 
obtain quadratic 
fit and approximate 
minimum AV # 


^ Is ^ 

180° transfer 
required 


Pass 

Number 


Figure II-3 OBJECT Program 
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Figure II-3 OBJECT Program 
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Figure II-4 COORDINATE SYSTEM AND TRAJECTORY SPECIFICATION 
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Appendix III (Unclassified) 

OPTIMAL-EXPLICIT GUIDANCE (OP-EX) FOR POWERED 
FLIGHT OUTSIDE THE ATMOSPHERE 

9 

This appendix describes the theory and implementation of OP-EX, 
an IBM-developed guidance algorithm for all phases of powered flight outside 
the atmosphere, including the exoatmospheric phase of ascent (i 0 e„ , that 
phase of ascent where the dynamic pressure is less than 30 p 0 s*f* ). As its 
name implies, OP-EX is both optimal and explicit,. That is, it accepts as 
inputs a given "present 11 state and propulsion performance models for the 
present and future stages of the rocket, and generates an optimal (fuel- 
minimal) trajectory v/hich satisfies explicitly stated constraints and final con¬ 
ditions, without dependence on off-board precomputations * OP-EX guidance 
is compatible with and is used by the Optimal-Explicit Flight Planner,, 

3. 1 SUMMARY OF CAPABILITIES AND ADVANTAGES 

OP-EX is a versatile guidance algorithm usable for exoatmospfreric 
ascent, orbit transfers, direct rendezvous, direct intercept, deboost - in fact, 
for all powered phases outside the atmosphere,, For each of these guidance re¬ 
quirements, the resulting fuel expenditure is minimal,. 

Also, for any given guidance phase, a variety of alternative terminal 
conditions can be easily specified,. For example, orbit insertions can be 
achieved with specification of (1) orbital plane and orbit orientation, size and 
shape in that plane, or (2) orbit inclination, latitude of perigee, size and shape, 
or (3) orbit inclination and latitude and longitude of orbit perigee* In fact, almost 
any reasonable combination of termination conditions imaginable can be easily 
implemented* All that is required is the specification of six "terminal-error" 
equations (according to a set of rules which will be derived) which completely 
define the mission type* These equations, when satisfied, determine the unique 
optimal trajectory* 

OP-EX can optimally perform any magnitude of plane change within 
the vehicle’s capabilities during ascent* It can also optimally perform any 
orbital maneuvers for plane change, orbit transfer, rendezvous, or intercept. 
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OP-EX has flexible provisions for propulsion performance descrip¬ 
tion (thrusts, stage masses, mass rates). Art unlimited number of stages can 
be specified and, in theory, any modelable function of time can be used to repre¬ 
sent engine mass flow and exhaust velocity. Even exotic schemes like the Saturn 
V PU (propellant-utilization) system, which involves a change of thrust level 
at a variable time after stage ignition, can be implemented satisfactorily. In 
fact, the optimal time at which to switch the level of thrust could be specified 
by OP-EX . 

The OP-EX algorithm is computationally efficient, to be practical 
for real-time guidance, yet is sufficiently accurate to be used for rapid predic¬ 
tion of AV requirements for flight-planning purposes. 

OP-EX provides predicted vehicle states at future staging times 
(and other critical times) which can be used in spent-booster impact prediction. 

OP-EX has growth capability for the optimal inclusion of additional 
constraints such as fixed final attitude, recovery ceiling, etc. Implementation 
would require little more than a modification of the control law; the routines 
for trajectory generation and for iterative adjustment of trajectory parameters 
would be essentially unchanged. 


3.2 


DERIVATION OF EQUATIONS FOR OPTIMAL TRAJECTORIES 


The equations of motion for rocket-powered flight outside the atmos¬ 


phere are 


r = v 


( 1 ) 


V = £(r,t) + a (t)u (t) 


( 2 ) 


where r^ is the position vector, v. is the velocity vector, j* is the acceleration due 
to gravity, a (t) is the magnitude of thrust acceleration, and u(t) is a unit vector 
giving the direction of thrust acceleration. For the one-burn case with non- 
throttleable engines, a (t) may be regarded as a given function of time, but in 
more complicated cases a propulsion-model must be included in the problem 
formulation. For an n-stage rocket, equations determining a (t) are 
a 


= T./M. 
3 J 




M. = -T./c*. + uj. \ 
J J J J / 


J 


It 2 p 


n 


(3) 


O < T < T 

j j max 


J 
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where M is vehicle mass, T is thrust, T max is full-throttle thrust, c* is the 
effective exhaust velocity, and cj is inert mass flow rate due to scheduled 
dumping of expendables and propellant leakage and evaporation. The effect of 
a; is, in general, very small, so to avoid needlessly complicating the dis¬ 
cussion, it will be neglected. (Should cu become significant, the theory can, 
of course, be expanded to include its effect. ) c*, T max , and M may change 
discontinuously at staging times. Also* c * and T max may be functions of the 
burning time since stage ignition. (For solid fuel rockets, this is usually the 
case.) The engines may be non-throttieable, non-restateable or both. How¬ 
ever, for formulating the optimization theory, it is convenient to begin by con¬ 
sidering perfectly throttleable engines. For cases in. which staging is triggered 
by fuel depletion, c* and T max may be regarded as given functions of the mass 
M. Also, (as a mathematical artifice) mass discontinuities at staging may be 
regarded as brief intervals with very small c* and correspondingly high values 
of M. 


Introducing a throttle variable S, 

a = Sa (t) 

max 

T (M)/M 
max 


max 
M 

O < S < 1 


-ST (M)/c*(M) 
max 


Equation (3) can be rewritten as 

(4) 

(5) 

( 6 ) 
(7) 


The optimization problem is to choose control policies u (t) and S (t) such that 
the resulting trajectory satisfies given final, conditions and extremizes some 
function of the final values of ^r, v;, M and t. Usually, but not always, the 
quantity to be extremized is the final mass M(t£). 

This is an optimization problem of the Mayer type. It can be ex¬ 
pressed in standard form by introducing the seven-component state vector 





LMJ 
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and the four-component vector 


v = 


Then the state equations can be expressed in the standard form 
x = _f (x, v, t) 

where 


(9) 

( 10 ) 


f 


"v 

£(r, t) + Su a (M) 
— — max 


(ID 


-ST (M)/c*(M) 

L max J 

The control constraints | u | =1 and O < S < 1 are independent of the state 

variables, so Pontriagin f s maximum principle is applicable. Introducing the 
adjoint vector 


£ 


ZZ 

X 


cr 


( 12 ) 


the Hamiltonian can be written as 

T 

H(x, £> t) = J£ (x> v, t) (13a) 

= 77‘V+X'g+Srx* ua -crT /c*l (13b) 
JL — — s L ~ max max J 

It is convenient to introduce the abbreviation K for the bracketed quantity in 

Equation (13), and G as the symbol for the gravity-gradient matrix. Then the 


standard adjoint equation 

£ = -dH/ dx (14) 

gives 

X = - 77 (15) 

£ = -G X (16) 

cr = -SdK/dM (17) 

The vector X is Derek Lawden T s ,f primer vector, " 

The maximum principle requires 

U = X/X (18) 
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and 


1 if K > O 


O if K < O 


An immediate consequence of Equations (6) and (17) is 
K = ( \ • u) a 


which simplifies to 

K =| X | a ' (21) 

1 1 max v 1 

when use is made of Equation (18). This relation enables K to be computed di¬ 
rectly, so the quantity cr becomes unnecessary and need not be computed. 

In many cases, the control of thrust magnitude is prescribed or 
trivial (full on till cutoff, zero thereafter) so only the optimization of steering 
need be considered. For such cases, the differential equations for an optimal 
trajectory reduce to 

£* = v = £(r_, t) +a(t)X/X (22) 

X = GX (23) 

For an inverse-square central field, the quantities g and G are given by 

£ = -(/j./r 3 )r_ (24) 

G = -(fi/r 3 ) [i - (3/r 2 )rr T ] (25) 

T 

where r_ is the transpose of the column vector jr. For multiburn cases, or 


cases in which the trajectory starts from a parking orbit with a freely chosen 
ignition time. Equations (22) through (25) must be supplemented by (19) and (21). 

Equation (21) shows that the quantities K and | \| always increase or 
decrease together, so local maxima or minima of K along an optimum trajectory 
are maxima or minima of |\|also. During coast phases, M and a max (M) are 
assumed constant, so K is linearly related to [ Xj . (Should an appreciable mass 
change occur during coast, due to scheduled dumping of expendables, propellant 
leakage or evaporation, or all three, the formulation could be expanded to include 
its effect.) Consequently, the magnitude of|Xj>must be the same at both ends of 
a co&st phase., .The powered phases occur at local maxima of | Mt)|. These 
maxima may occur at the beginning or end of the trajectory (exterior maxima) 
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or at intermediate points (interior maxima). For interior maxima (i.e., burns 
whose beginning and end times are not determined by any consideration other 
than optimality), the conditions for optimal thrust control require the integral 

J = / I X 

burn 

to vanish. This requirement does not hold for a burn whose beginning time or 
end time is constrained. 

The optimality conditions given above reduce the problem of tra¬ 
jectory optimization to a two-point boundary value problem (TPBVP). The ini¬ 
tial conditions generally consist of given values of r_(t Q ), v(t Q ), M(t Q ), and t Q . 

For cases in which the final mass is to be maximized, the final conditions con¬ 
sist of k < 6 mission conditions {prescribed relations among the final values of 
v(t£), and t£ plus 6-k transver sality conditions which come from optimal 
control theory. The transver sality conditions can be analytically derived from 
the mission conditions by the requirement that 

k ‘ - h ' = G t - t f (2?) 

must hold for every pair of infinitesimal variations _Sr, Sv consistent with the 
mission conditions and at every point on the coast trajectory. 

3.3 NUMERICAL SOLUTION OF THE OPTIMAL TRAJECTORY PROBLEM 

As derived in the previous section of this appendix, the differential 
equations defining the optimal trajectory of a rocket powered vehicle are 

£ = v = g(r) + a(t) u (28a) 

X = GX (28b) 

u = x/x (28c) 

A guidance scheme also based on Equations (28), but differing from OP-EX in 
computational techniques, has been described by Brown and Johnson (Reference 
1). Numerical integrating of these equations from some initial time t Q to an 
estimated final time t f , starting from the initial state rjt 0 ) and v(t Q ) and esti¬ 
mated initial costate X(t Q ) and X (t Q ), determines the final state, r/tf) and vftf), 
and the final sqstate variables X (tf) and X (t£). Seven error quantities, derived 
from k specified final state conditions and 6-k transver sality conditions and a 


max 


dt 


( 26 ) 
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desired normalization of X(t Q ), must be brought to null by iterative adjustment 
of the seven parameters X(t Q ), X (t Q ) and t f . To this end a Newton-Raphson 
iteration procedure will be used. 

Since the same trajectory integration algorithm is to be used for 
both flight planning of all powered phases and for operational guidance during 
the powered phases, the integration must be both fast and accurate. An obvious 
method of increasing the speed of solution is to increase the step-size used in 
numerical integration of Equations (28). Since the guidance is rr closed loop rr and 
can be made accurate near cut-off at small cost, this does net cause guidance 
errors. However, it has three possible adverse and unavoidable effects: 

(1) The steering policy becomes non-optimal, incurring a 
fuel penalty. 

(2) The initial estimate of t^ becomes inaccurate, which distorts 
flight planning. 

(3) Truncation errors in the computation of sensitivity coefficients 
(which are required for iterative correction of the estimated 
parameters X (t Q ), X (t Q ) and tf) cause them to differ from the 
actual sensitivities of the numerically integrated Equations (28). 
This degrades the convergence of the iterative solution, and 

in extreme instances, destroys convergence completely. 

However, it is known that present day explicit guidance schemes, 
which are essentially one-step integrations of the trajectory equations improved 
by closed formulae for certain thrust integrals, have small performance penal¬ 
ties. It may therefore be expected that an "explicit-predictor tf integration 
scheme which resembles the application of explicit .guidance to each integration 
step will have negligible performance penalty for very large integration steps. 
This, in fact, is the case. It has been demonstrated that one integration step 
per stage of the Titan IIIC booster produces a AV performance penalty of only 
0.1 ft/sec. 

Errors in initial estimates of t£, made when the time remaining un¬ 
til thrust cut-off is large, are quite significant for explicit guidance schebaes, ' 
but have been found to be acceptably small when the explicit-predictor integration 
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scheme is used with large steps. It has been repeatedly demonstrated, using 
the Titan IIIC booster, that integration steps o*f 100 seconds produced an error 
in predicted cut-off time of less than 0.3 seconds in cases where the time re¬ 
maining until cut-off was greater than 700 seconds. Even larger time steps 
are believed useable for mission planning. 

The remaining factor governing useable step-size is the effect of 
truncation error on the convergence of the iterative solution. This effect dis¬ 
appears if the sensitivity coefficients are computed by finite differences. This 
is done by first generating an n unperturbed M trajectory and evaluating the error 
quantities. Then fr perturbed Tr trajectories are generated by changing components 
of the initial adjoint vector one at a time by a small specified variations. Sub¬ 
tracting the error quantities of each perturbed trajectory from the corresponding 
error quantities of the unperturbed trajectory, and dividing by the change in 
the initial adjoint component, gives one column of the sensitivity matrix. 

Computing sensitivities by finite difference is similar to performing 
experiments on the numerical trajectory solution. Regardless of truncation 
errors, approximations in the equations, etc., the sensitivity coefficients gen¬ 
erated by finite differences will give an accurate prediction (except for round¬ 
off and nonlinearity effects) of the way in which the numerical trajectory solu¬ 
tion will change if its parameters are changed. 

Experience with the ascent case has shown that there is no difficulty 
in choosing parameter changes large enough to avoid trouble from roundoff 
errors, but small enough for linearity to hold. In trial experiments, the sizes 
of the parameter increments were varied over a range of more than three or¬ 
ders of magnitude without appreciable degradation of convergence. 

3.4 EXPLICIT-PREDICTOR INTEGRATION ALGORITHM 

The "explicit-predictor 11 integration algorithm is designed to take 
advantage of closed formulae for integrals involving thrust acceleration. The 
formulae presented here assume that rocket engine mass flow rate M and effec¬ 
tive exhaust velocity c* are constant. This should not be construed as a limita ¬ 
tion of the technique . In theory, at least, any cpriodelable function of burning 
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time can be used for either M or c*. For example, a n propellant-utilization" 
model in which the level of thrust is altered sometime after engine ignition 
could be implemented within the framework of the algorithm. Also there is 
no limit imposed on the number of rocket stages. 

The equations of the integration algorithm will be shown for the jth 
rocket stage and the nth time step. Since engine parameters and mass are 
generally discontinuous at staging, the integration time step is selected so that 
no more than one stage is ever included in any integral evaluation. In other 
words 


At = Min (At , t..- t ) 
n max j n 


n+1 


s t + At 

n n 


(29) 


where t- is time at the end of stage j, t is present time and At ^ is the 
j ° J n max 

largest time step allowed. 


The thrust acceleration term a(t) u is separated into two parts, that 
which defines its magnitude a(t) and'that which defines its direction u. Over 
each integration step, u is approximated by a second-order vector function of 


time, 


u(t) = u (t ) + a (t-t ) + P (t-t ) t < t < t 

~~ m “n m —n m n 


n+1 


(30) 


where 


t = t + At /2 

m n n 


As will be shown subsequently, it is possible to predict u(t ) an d n(t^_^). Given 
the three values of u, the vector coefficients are 


a = F u(t ) - u(t )] / At 

— n L n+1 — n J n 

p = 2 [u(t ) + u(t ) - 2 u(t )] / At 2 . 

—n n+1 — n — m J n 


(31) 



Am-io 


from t to 
n 


where 


M . is the 
°J 


To complete the single and double integrals of thrust acceleration 
t n + 1 * ^ i- s necessary to evaluate six integrals: 
t 


I. 


n+1 


C' 1 '. 


n 


T. + t. - t 
J J-l 


n+1 

c*. (t-t ) 

_J m 
r. +~t. , - t 
J J-l 


dt 


dt 


n 


n+1 

1° c*. (t-t Y 
I _j _ m 
J t. + t - i 

& 3 J-l 


n 


n+1 


dt 



dsdt 


(32) 


n+1 



c*. (s-t ) 
_ 3 m 

T. + t. . - S 

J J-l 


dsdt 


t t 
n n 


n+1 t 

0 2 
c*. (s-t ) 

_J._EL_ 

r. + t. - s 
t t J J-l 

n n 



dsdt 


r = -M ./M. 
j oj J 

mass at the igmtiori.of stage j. ^ 
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stant 


For the assumptions made above, namely that c* and M are con- 

J j 


c* log (1 + p ) 
J e 


where 


c* At - ( T. + t - t ) I 
J n J J-l n+1 1 


~ l A - & ~ t 4.1 > *t 

4 m n+1 1 


( T. + t. , - t ) (I - C* At ) 

J J-l m 2 j n 

+ c* ( r. + t. - t ) 2 (2 p + p 2 )/ 2 
J j j-l n +1 r r 


( t. + t. - t ) I - c* At 72 

j j-l m 4 j n 

( r. + t. - t ) (I,. - c* At ^ /2) 

j j-l m o j n 

" C *j [ ( T j + Vl " t n +l )3 (3 ^ + 3 P + ^ 3)/6 

- ( T. + t - t ) 2 At / 21 

j J-l n n J 


At / ( r. + t - t ) 
n j j-l n+1 


(33) 


Having constructed the closed formulae for the first and second 
integrals of thrust, the total "explicit-predictor" algorithm can now be shown. 
The contributions to position and velocity due to thrust acceleration, Ay f and 
A v 1 , and due to gravitational acceleration, Ay" and Av",will be isolated as a 
hedge against requiring extended precision integral accumulation. The primer 
vector can now be predicted at t m and by using the second order equations; 

X = \ + \ At /2 + \ At 2 /8 

m n n n n n 

U = x / x 

“~m *~m m 


(34) 


X . = X + X At + X At /2 

“"n+1 ~n —n n —n n 


u = x /X 
•Vl-l -n+I ' n + 1 
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After evaluating a and ft according to Equations (31) and 1^, * * *, 1^ according 
to Equations^(33), the thrust acceleration contributions to position and velocity 
are: 

+ 1. s„ + h A. 


Ar 1 , = 

Ar' 

+ Av 'At 

— n+1 

~ n 

— n j 

Av 1 = 

Av' 

+ I u H 

- n+1 

n 

1-m 


A 

U 


(35) 


-n 


The gravitational acceleration contributions to position are: 


A r. ,f 


n+1 


Ar !r + Av rr At + g At 2 /2 + G v At ^/6 
~~ n n n ni n n ~n n 


(36) 


Total position is 

r T = jr + v (t . - t ) + Ar_ ! . + Ar f? . . 

^i+I o o n+1 o nri ““ n+1 (37) 

Gravitational acceleration and the gravity gradient matrix can now be evaluated 

•* 

using the predicted value of r_ n _j_p then X n+ -^ can be evaluated 


4.+1 = - ( ^ /r n+l )x n+l 
G n+I =-W r n+1 3 >[l-(3/ 


2 T' 

'n+1 -n+1 —n+1 . 


(38) 




^n+1 ~tl+1 


Finally, velocity and primer vector rate are 

X = X + ( X + X ) At /2 

~“n+l “~n ~n “"n+1 n 

Av” +1 - Av n + (c + g ) At /2 

n+1 n =n+l ~n n 


(39) 


X. 


n+1 


v + Av n + Av 1 ,, 
o — n+i — n+I 


Equations (34) to (39) are the explicit-predictor integration technique. As the 
problem is integrated, the index j must refer to the lf current f - stage. It is worth 
noting that no restrictions are placed on the time at which the problem is begun. 
In this way the algorithm is applicable to the real-time guidance requirement. 


3.5 


ITERATIVE DETERMINATION OF ESTIMATED PARAMETERS 


As pointed out previously, seven terror quantities must be brought 
to null by iteratively adjusting the estimated parameters X (t Q ), X (t ) and t£.' 
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Also, since the normalization of the primer vector at t Q is irrelevant, by an 
expeditious choice of coordinates a "preferred coordinate" can be established 
along which the component of X can be left invariant. What this means is that 
the "preferred axis" should not be chosen orthogonal or nearly orthogonal to 
the expected direction of thrust acceleration. In general, this is an easy re¬ 
quirement to meet and is worthwhile since it reduces iteratively corrected 
parameters from seven to six. 

It is now necessary to establish the six-by-six matrix of sensitivities, 
the so called Jacobian matrix, which relates the variation of estimated param¬ 
eters to the variation in error quantities. The error quantities used are unique 
to, and indeed define, the type of trajectory being generated. Error quantities 
for several trajectory types will be presented in a subsequent section. 

The Jacobian matrix will be computed, as was stated previously, by 
employing finite difference techniques. Only five "perturbed" and one "unper¬ 
turbed" trajectories must be generated. (A*t this point, one can see the obvious 
advantage of an extremely rapid trajectory generation technique.) Each per¬ 
turbed trajectory yields one column of the Jacobian matrix. The sixth column, 
which represents the changes induced in the error quantities by a unit charge 
of the final time, can be evaluated by extrapolating the final state of the unper¬ 
turbed trajectory over some small time increment and then reevaluating the 
error quantities. 


If the error quantities for the six perturbed trajectories are desig¬ 
nated by column vectors . . . , of dimension 6, and the unperturbed tra¬ 

jectory error quantities by c Q , the six-by-six Jacobian matrix is 


= \(£ - € )/AX ( c - eJ/AX (c - O/AX , ( <= - O/AX , 
”~o ““1 , o 2 , o 3 o 4 


€ - O/AX .. ( € - € )/ At 1 

O 5 O 6 f J 


J 
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The Newton - Raphson iteration equation in vector form is 


-1 

^ • 

1—> 

_1 



*•2 


^1 

*3 

- 

*3 

*2 


^2 

x 3 


^3 

J- 

Hn 

1_ 


_ t f_ 


( 41 ) 


where the numbered subscripts on X and X refer to the first, second, and 

third components. Note that X^ is arbitrarily chosen to be along the "preferred" 

coordinate. 

In Equation (41), F is a convergence progress control parameter 
which is adjusted according to a "performance indicator". The performance 
indicator is defined as 


P 


6 

2 w. 

i=l 1 


2 

€ . 
Ol 


(42) 


which is the weighted sum square of the unperturbed error quantities. The 
weighting factors W., are required because the error quantities are not (in gen¬ 
eral) of the same dimensions. 

It can be argued that if the Jacobian matrix is correct and non sin¬ 
gular, progress must be made by Equation (41) in reducing the error quantities 
if a sufficiently small T is used. On the other hand, if T is too small, the re¬ 
duction process will be unnecessarily sluggish. The rule used to adaptively ad¬ 
just r is as follows: if the performance indicator on the present iteration is 
less than it was on the last iteration, T is set 

T —Min (2r, 1) 

and the evaluation of a new Jacobian proceeds. If the performance indicator is 
greater than it was on the last, iteration, T is set 

r — r /2 
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3.6 SPECIFICATION OF ERROR QUANTITIES 

As shown previously, it is necessary to specify six error quantities 

9 

which uniquely define a mission type and, when brought to zero, satisfy the 
mission objectives. Of these six, k are specified final state conditions (or func¬ 
tions of final state conditions) and 6~k are transversality conditions which in¬ 
sure an optimal trajectory. 

Three basic options most used by QRGT will be described below. 

The other orbit injection options described in the summary of capabilities can 
be derived in a straightforward manner but are not presented here. 


1) Injection into Orbit, with Specification of Orbit Size, Shape, and 

Orientation. 

For injection into orbit with prescribed orbit size, shape and orien¬ 
tation, the only degrees of freedom remaining are the time and location at which 
orbit insertion is to occur. Let us assume that the desired orbit is specified by 
its radius and velocity vectors at perigee, rp and v . With this specification, 
the required velocity vector at any location in the orbit is given by 



where Q is the true anomaly; The semilatus rectum, s, is given by 


2 2 / 

s = v r /u i. 

P P/^ 

The magnitude of the required radius at any point in the orbit is 


(46) 


_ 5 

r R (1 4- e cos Q ) 

where the orbit eccentricity is given by 


1/2 


(47) 


e 


(s/r - 1) 
P 


(48) 
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The unit vector normal to the orbital plane is 



P P 


(49) 


As shown above, the required velocity vector and required radius 
are functions of the true anomaly. The true anomaly of the generated trajectory 
is given by 


e 


-i 

- Tan 


r.(tJ * v /v 

{ _P P , 

r.ftj * r. fr 
f p p 


Five of the six required error quantities are 


(50) 


€ n5 

€ nG 


= x R (9) - v(t f ) 


e n2 


= r R (9) - r(t f ) 


(51) 


£ n3 ' "a ’ *<*£> 

These five quantities are sufficient to satisfy all mission specification. The 
sixth quantity gives the deviation from satisfaction of a transversality con¬ 
dition, which is 


A (tf) *£(tf) -^(tf) *v(t f ) = 0 (52) 

and expresses the fact that phase-in-orbit is not specified. This is the well- 
known transver sality condition for n time-free ,r cases. 


2) Direct Ascent to Rendezvous, with Rendezvous Time Specified 

Direct-ascent-to-rendezvous means that, after thrust is initially 
terminated, only a velocity matching burn at the specified rendezvous radius is 
required. For this case, the required velocity, v-^ is determined by means of 
Lambert f s theorem. One statement of Lambert r s theorem is that, if two posi¬ 
tion vectors, the time required to traverse between them, and whether the true 
anomaly to be traversed is greater than or less than 180° are specified, the 
orbit of traverse is unique and deterministic. The equations used to determine 
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can be found in Reference 2 and are not repeated here. The remaining 
three error quantities must be specified by transver sality conditions. 

To determine the applicable transver sality conditions let us first 
examine the equation of perturbed motions obtained from Equation (28a). 

8V = Sv = G Sr. + 8a (53) 

To obtain Equation (53), it is assumed that | 8 r. | is small enough to justify ig¬ 
noring terms of the order of |8rJ^. 

Recalling Equation (28b), we see that on a coast trajectory, where 
Sa = 0, Equation (28b) is identical in form to Equation (53). Hence, the well 
known state transition matrix solution for Equation (53) must also be a solution 
for Equation (28b). Symbolically, 


"X(t: r r 

X(t r ) 


$ (t , t ) 

r f 


_ x(t f r 

X(t f ) 


(54) 


where t r is the time of rendezvous. If \(t r ) is assumed to be zero (i.e., a fixed 
attitude orbit matching burn) the required primer vector rate at thrust cut-off 
is given by 


where 


X R (t £ ) = ® 




-1 

X(t ) - 

12 

l r 

~<$ 11 

I 12“ 

21 

! $ 22_ 



(55) 


All that remains is to determine \(t r ). Recalling Equation (21), we 
see that for an optimal trajectory the magnitudes of the primer vector at the 
first thrust termination and at the beginning of the velocity matching burn must 
be equal. The direction of the primer vector at the beginning of the velocity 
matching burn must be approximately in the direction of the velocity-to-be- 
gained vector. Symbolically, 

x (t) = x(tj i^rr^r ■ (56) 

r I L 1 r r J 
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where v.rp (t r ) is the target velocity ^ndv^t^.) is the interceptors velocity at ren¬ 
dezvous. Hence, the three additional error quantities are 

£ nl 
€ n2 
/n3 

Reference 3 describes a very simple and elegant way of computing 
the required state transition matrix which is valid for all orbits except the de¬ 
generate rectilinear orbit. 

3) Direct Ascent to Intercept--Intercept Time Specified 

Direct-ascent~to~intercept is very similar to direct-ascent-to- 

rendezvous just described. The principal difference is that no velocity matching 

burn is performed at intercept. The required velocity for intercept v\, is com- 

.K 

puted in the same manner as for rendezvous but the transversality conditions 
are not. 

Recalling Equation (2 7), we have 

X * Sr = X • Sv t ^ t£ (58) 

This equation must hold at all points on the coast arc. If we choose the point at 
intercept, the perturbation in position, Sr^, must, by definition, be zero. Since 
there is no constraint on velocity at intercept, the only way Equation (58) can be 
satisfied is if X at intercept is also zero. 

Equation (57) can now be rewritten for the intercept case and is 

-1 

= - $ <£ X(t ) - X(t ) (59) 

12 11 

3.7 USE OF OP-EX FOR ORBITAL MANEUVERS 

OP-EX has several attributes which make it ideally suitable for 
guiding orbit maneuvers as well as ascent. Principal among these is its ability 
to optimally perform large plane changes. In fact, the direct-ascent-to-rendezvous 


€ nl 

e n2 

€ n3 


-1 


= $ 


12 


X(t f ) 


in (t ) - V(t ) 

1 r r 


~ $ X(t)UX(t) (57) 

n i I 
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option can be directly applied (except for minor changes in initialization) to all 
intermediate burns and their subsequent coast arcs. As has been shown pre¬ 
viously, additional guidance options can be easily added without sacrificing the 
optimal performance of OP-EX. 

In keeping with the QRGT mission planning policy, the time of ini¬ 
tiation of each thrusting maneuver and the location and time of the subsequent 
thrusting maneuver are determined by the mission planner. Since the planner 
assumes impulsive changes in velocity, the actual time of engine ignition during 
active guidance will be modified to make the nominal centroid time of the burn 
correspond with the impulse time specified by the planner. 


3. 8 EFFECT ON OP-EX OF MEASUREMENT OF ENGINE PERFORM¬ 

ANCE DUR IN G ACTIV E GUIDANCE 

OP-EX guidance relies on propulsion performance descriptions for 
the various stages to extrapolate the trajectory forward to thrust termination. 

In active guidance, variation of engine performance from that specified will 
cause the extrapolated trajectory to be in error. However, because the guidance 
is ff closed-loop ,f by design, (that is, the trajectory is recursively extrapolated 
every ten seconds during early stages and more frequently in the final stage) 
such performance variations do not seriously effect the satisfaction of guidance 
objectives, but the fuel optimizing ability of OP-EX is handicapped (although 
certainly not disastrously) by the variations, 

A ,r mass-rate n prediction scheme, which will be derived subse¬ 
quently, was implemented for study purposes and was found to enhance the fuel 
efficiency of OP-EX guidance in the presence off-nominal engine performance. 
The amount of AV saved, in every case, was a few ft/sec. 

The ft mass-rate ir predictor is formulated as follows: the theoret¬ 
ical AV gained in stage j between the k and (k+1) predictions of OP-EX is (if 
c* is constant) 







AIII-21 


where M , is initial mass, c*. is effective exhaust velocity, t. . is time of 
J J J ~ 

ignition (a measured quantity) and M. is estimated mass rate of the jth stag 

J 

The measured AV 1 gained is approximately 
AV'=£(t k )- [w(t k+1 ) - w(t k >] 


( 61 ) 


where u(t) is the direction of thrust acceleration and w(t) is the velocity-meter 
output vector at time t. Predicted average mass rate is then 




M. = 
J 



- 1 - exp AV'/c* jj 


-t. 1 - T t - t 1 exp f- AV'/c*. \ 
j-lj L k j-lj r \ • J J 


( 62 ) 


It can be seen from Equation (62) that errors in initial mass are 
absorbed in the estimated mass rate. 

If every last ounce of fuel is to be saved, perhaps a multistate filter 
of the Kallman type could be used to estimate M -, M. and c*-. 

J J J 
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GROWTH CAPABILITY FOR TREATMENT OF CONSTRAINTS 


For some vehicles and missions, the flight guidance must satisfy 
constraints other than the usual constraints on final position and velocity. Ex¬ 
amples are (1) recovery-ceiling constraints and (2) requirements for a pre¬ 
scribed attitude at cutoff, together with a low angular rate. The structure of 
the OP-EX guidance algorithm facilitates the introduction of modifications to 
make it capable of handling such constraints. The principal changes required 
are (1) modification of the form of the control policy, in a way that introduces 
new trajectory parameters, and (2) providing equations for computing new 
error quantities which measure constraint violation. The explicit-predictor 
algorithm of OP-EX can be used, without change, to generate trial trajectories 
by which the old and new error quantities can be evaluated for any given values 
of the old and new trajectory parameters. 

Before active guidance begins, the on-board flight planner (using 
part of OP-EX as a subroutine) will generate a feasible, optimized nominal tra¬ 
jectory and so provide an initial set of values for the trajectory parameters. 

Some of these trajectory parameters will remain constant during active guidance, 
others will be iteratively adjusted by OP-EX. For this purpose, the matrix of 
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sensitivity coefficients must be enlarged by additional columns corresponding 
to new error quantities, and additional rows.for new trajectory parameters. 
Since the matrix elements are evaluated by finite-difference methods, this en¬ 
largement is straightforward. 

Appropriate forms for the control-policy modifications can be de¬ 
rived from optimal control theory. It is not necessary to implement the true 
optimal solutions (which would require, in general, the computation of new 
adjoint variables) since it is generally not difficult to find simple approxima¬ 
tions which give satisfactory performance. 

For example, for the recovery-ceiling constraint, a preliminary 
analysis indicates that a satisfactory scheme may be realizable as follows: 
the control policy is modified by introducing discontinuous changes of \ and X 
at a scheduled future time t g which is prior to the predicted cutoff but after the 
latest time at which the recovery-ceiling constraint might be active. These 
changes are of ^he form 


AX = kb^ 
AX = k^b £ 


(63) 


where k is a parameter to be iteratively adjusted by OP-EX. The vectors b^ 
and b 7 are generated by the flight planner, together with the time t g , and are 
not changed during active guidance. When the danger of violating the recovery¬ 
ceiling has passed (which will occur before time t g ) the parameter k is set to 
zero, restoring the normal control policy. At times immediately after this 
switch-over, the rate of change of commanded thrust direction must be limited 
to prevent undesirable transient behavior. 

The problem of satisfying attitude constraints at cutoff involves both 
guidance and control. Optimal control theory indicates that the optimal trajec¬ 
tories consist of two arcs. On the first arc, which covers most of the trajectory, 
the present OP-EX control policy is optimal, and constraints on angular rates 
or angular accelerations have no direct influence. On the second arc, the atti- 
fcude control system is working at maximum effort to achiey§ the desired final 


attitude and attitude rate. In an actual implementation, the maximum-effort 
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phase must be followed by a precision-control phase that tries to realize 
accurate end conditions. A reasonable approach to the over-all problem is to 
first consider the attitude-change maneuver as a separate optimization prob¬ 
lem (intimately related to the autopilot design) and devise a quasi-optimal 
policy (with margins for performance uncertainties) for achieving a given final 
orientation with a low final angular ra.te. This gives a formula for the maneu¬ 
ver-time required, as a function of augle-to-be-turned-through, which is valid 
when the initial angular rate is low. This formula (and approximate forms for 
the angular rate profile during the maneuver) can be used in the explicit- 
predictive algorithm for generating predicted trajectories. When the attitude- 
change maneuver begins, control of cutoff position will probably have to be re¬ 
laxed. However, accurate control of both attitude and velocity is possible by a 
policy which rapidly nulls the components of velocity-to-be-gained which are 
normal to the desired final roll-axis direction, and controls the parallel com¬ 
ponent by cutoff time. 

3.10 SIMULATION RESULTS 

The OP-EX guidance algorithm has been tested by simulation, using 
the IBM-developed GISMO simulation program, with the following guidance 
options: 

1) Direct ascent to intercept; initial ascent plane coincident 
with target plane 

2) Direct ascent to intercept; target plane rotated 40° from 
initial ascent plane 

3) Direct ascent to rendezvous; initial ascent plane coincident 
with target plane 

4) Direct ascent to rendezvous; initial ascent plane rotated 40° 
from target plane 

5) Insertion into circular orbit; no dog-leg during ascent 

6) Insertion into circular orbit; 40° dog-leg during ascent 

7) Insertion into noncircular orbit (eccentricity 0.047) with 
no dog-leg during ascent 
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8) Insertion into noncircular orbit (eccentricity 0. 047) with 
40° dog-leg during ascent 

• 9) Direct rendezvous from orbit; 40° plane change 

10) Orbit insertion from orbit; 40° plane change 0 
In all cases, the Titan IIIC booster was used. Three payload options 
were used for the various missions. A spherical Earth Model in GISMO was 
used for all tests since gravity oblatness correction equations were not complete 
at the time of testing. 

In every case the guidance algorithm performed with aplomb. No 
evidence of non convergence was indicated. The number of iterations for initial 
convergence ranged from one for option 4) to nine for option 1). This difference is 
probably due to the fact that all ascent cases were initialised with the same 
standard set of initial parameter values, and this set is closer to the requirements 
for injection into a circular orbit than to those for intercept. For all orbit-to- 
orbit burns, the primer vector is initially aligned along the estimated velocity-to- 
be -gained vectox furnished by the mission planner. It is possible that this rule 
(or an adaptation of it) should be followed for direct-ascent intercept cases also. 

In all intercept or rendezvous cases, the error in position at the 
specified time of intercept was less than 20 feet. In all orbit insertions, the 
error in radial and out of plane position components was less than 20 feet, and the 
error in velocity was less than 0.05 ft/sec. These errors were, of course, for 
nominal engine performance and probably represent, to a high degree, the numeri¬ 
cal resolution of the 7094 computer. 

3.11 OP-EX MATH FLOW DIAGRAMS 

Math flow diagrams of OP-EX guidance are presented in Figure III-l. 
Table HI-1 gives math flow symbol definitions. 




Set Newton** 
Raphson 
Iteration Gain 
K = 1 


OP-EX GUIDANCE ALGORITHM 
Major Loop Entry 


► Fir stV ^ Yes 
Entry^ 


Go initialize 
for 

present mode 


AIII-25 

"ftTII 


yr HasN. 
'Stage Change' 
Ss sPccur r e d S' 


. .Is ' 
M(i) = 0 


Compute Velocity 
Meter Increment 

" J Q * 


mL ; 



















Go Generate' 
Trajectory 
and Error 
Quantitie s > 


P4-1 


o 

0 

mpute 

Column k] 

of 

Jacobian ! 

J," 

= \e ~ 

el/AX 

-«*k 

L-o 

-J o 

X' 

(k)-X* 

(k)-AX 

o 

c 

) O 




















Iteration Converged 

Enter t in Interrupt 
co 

Queue for Engine 
Termination 


Res 

et Time~of- 

Cut 

-Off 

; t 

•*- t (m- 1} + 0. 1 

CO 

f 
















Ain-2 8 











AIII-29 


At 1 <C) 


Return to 
Calling 
Routine ^ 


Compute Thrust Integrals 

X- - M 1 (m)/M l (m) 

Q 

p~ At 1 / Tt+ ~ t ! J 

I x = C(m) log^ (1 +/>) 

1^ - C(m) A t ! ~ jV + t^(m-l) -» Ij 


1 = ~ I ~ ft - t*T I 

2 4 [ m j 1 


t^(m~l) ** t^j ^ ~ C(m) At*j 
+C(m) t^(m«l) - t r j j^2p 4 /T'j/2 

I = r+ t r (rn~l) ~ t T ~ C(rn) At 1 /2 

5 L f mj 

I = T+ t (m~l } - t [l ~ C(m) At j2 /2j 

6 L f ml L 5 J 

- C(m) ^Jr + t f (m~i) ~ t*J 3 2 + 3/>j / 

- jr + t^(m^l) ~ t 1 4A t^J At ,2 / 2 } 


X = X + X* A t x /z + X, At 1 /8 


♦m ■♦co -♦co 

u = X / X 

—►m rn 


X — X -t X At* -t X* At 1 /2 

-♦CO “♦CO -♦co .-♦l 

u = X /X 

-♦2 -♦co co 


r = r + vAt* 

rj-i B r^i m - nn 


Figure III-1 OP- EX Guidance Algorithr 


P6-1 



















Set ‘‘Initial-Phase 1 * and 
'‘Last-Phase 1 * Indices 


i = 1 


Set up n Engine 

-Li 

Descriptions 

M, M , t, 

• o d 


Return to 
Calling Routine 


Figure III-1 OP-EX Guidance Algorithm 


Sheet 7 of 9 









* Error Quantity Generation Subroutine AIII-32 


Option 1 



Figure III-1 OP-EX Guidance Algorithm Sheet 8 of 9 













Ain-34 


C(m) 

e 


1.1 s 1.2 


G 

o 

G 


i 


I I* • * * sI 6 

j 

J 

J 

♦n 

K 


M (m) 
o 

M' (m) 
o 

M(m) 
M ! (m) 



m 

n 

A 

n 


Table III-l 

OP-EX GUIDANCE MATH FLOW 
SYMBOL TABLE 

Effective exhaust velocity of stage m 

Eccentricity of specified orbit (option dependent) 

Gravitational acceleration at beginning and end of trajectory 
integration step 

Gravity gradient matrix at t 

Gravity gradient matrix during integration 

Initial stage index 

Thrust integrals used during integration step 
Newton-Raphson iteration count 
Jacobian matrix of sensitivity coefficients 
nth column of Jacobian matrix 

Newton-Raphson convergence control coefficient 
Nominal initial mass of stage m 
Estimated present mass of stage m 
Nominal mass flow rate of stage m 
Estimated mass flow rate of stage m 
Nominal fuel mass in present stage 
Present stage index 
Final stage index 

Unit vector normal to specified orbital plane (option 
dependent) 
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Table HI-1*- OP-EX Guidance Math Flow Symbol 
Table (continued) 

Unit vector along desired body-fixed roll axis 
Present value of performance indicator 
Last value of performance indicator 

Required value of performance indicator for successful 
convergence 

Unit vector along desired body-fixed pitch axis 
Position vector at t 

o 

Gravitational acceleration contribution to position vector 

Thrust acceleration contribution to position vector 
Position vector at t 

co 

Position vector at intercept or rendezvous (option dependent) 
Specified perigee position vector (option dependent) 

Semilatus rectum of specified orbit (option dependent) 
Present time 

Time at start of trajectory generation 
Present predicted time of thrust termination 
Last predicted time of thrust termination 
Nominal burning time of stage m 
Estimated remaining burning time of stage m 
Estimated time at end of stage m 
Time of last guidance iteration 

Time of intercept or rendezvous (option dependent) 
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Table III- 1. OP-EX Guidance Math Flow Symbol 
Table (continued) 


t 

m 

t 

po 

t J 
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AAA 
U, , U , U 

~*m -»2 
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■**o 



v 

-rm 

-*mL 

W k 

«. a 

8t 

CO 

sx 

—~o 


Time at midpoint of trajectory integration step 
Time at beginning of present stage 
Running time during trajectory integration 
Steering vector at t 

Steering vectors at beginning, midpoint and end of trajectory 
integration step 

Velocity vector at t 

Gravitational acceleration contribution to velocity vector 
Thrust acceleration contribution to velocity vector 
Velocity vector at t 

co 

Target velocity at rendenzvous (option dependent) 

Interceptor velocity at intercept (option dependent) 

Specified perigee velocity vector (option dependent) 

Required velocity (option dependent) 

Present velocity meter output 
Last velocity meter output 

Weighting factors for performance indicator computation 

Coefficients of parabolic approximation of steering vector 
over one trajectory integration step 

Correction to estimated time of thrust termination 

Correction to estimated initial primer vector 
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Table III- 1. OP-EX Guidance Math Flow Symbol 
Table (continued) 

Trajectory integration time step 
Maximum trajectory integration time step 
Perturbation value for t 

co 

Velocity meter increment added since last guidance iteration 

Perturbation value for components of X 

Error quantities for unperturbed trajectory 

Error quantities for perturbed trajectories 

True anomaly in specified orbital plane (option dependent) 

Estimated primer vector at t 

Primer vector at t 

co 

Last estimated primer vector at t 

Primer vector at midpoint of trajectory integration step 
Gravitational coefficient of spherical Earth model 
n defined by equation 11 
,f defined by equation" 

State transition matrix relating state at t c to state at t. 

(option dependent) 

Three-by-three partition of state transaction matrix 
(option dependent) 
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Appendix IV 

.DIRECT SEARCH OPTIMIZATION PROGRAM "DSOP" 

1. INTRODUCTION 

A FORTRAN IV program named DSOP (Direct Search Optimization 
Program) has been developed for the experimental study of optimization algor¬ 
ithms and their application to mission planning in the QRGT study. DSOP finds 
a local minimum of any given function of a set of independent variables, sub¬ 
ject to a given set of constraints. Explicit constraints (bounds on the independent 
variables) are handled directly. More general constraints are treated by the 
penalty-function method. 

As the name implies, the optimization algorithms of DSOP are of 
the ,! direct search 17 type; that is, they make no use of analytical derivatives. 

A test problem (or any problem to which DSOP is applied) is therefore com¬ 
pletely defined by equations for evaluating the function to be minimized and the 
constraint functions (if any). No equations for partial derivatives of these 
functions are required. This characteristic is highly desirable in an optimizer 
designed for use on-board a spacecraft, since it greatly reduces the number of 
routines that must be provided if there are a large number of different optimi¬ 
zation problems which may require solution. Also, it makes the optimizer 
applicable to problems where analytic derivatives would be difficult or expen¬ 
sive to generate, and improves flexibility by simplifying the specification of 
new problems for optimization and the modification of existing ones to accom¬ 
modate changes in mission objectives and/or constraints, etc. 

DSOP consists of a main program and ten subroutines. This struc¬ 
ture facilitates experimentation since the interconnection of the subroutines can 
easily be changed. The structure of MAIN is described in Paragraph 2. 

The two principal search methods used are PMS (Pattern Move 
Search) and VA04A (Powell’s method). PMS is a modification of the direct 
search technique developed by Hook and Jeeves (Ref. I ). VA04A is described 
in Powell’s 1964 paper (Ref. 2). PMS, VA 04 A, and the other search methods 
used are described in Paragraphs 3 and 4 of this Appendix. 
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Results of this optimization experiments on test functions are pre¬ 
sented in Paragraph 5. For the experiments, without constraints, PowelPs 
method (VA04A) gave the best performance. Experiments using penalty func¬ 
tions to enforce constraints were successful with PMS, but relatively unsuc¬ 
cessful with VA04A # The reasons for this are not fully understood at present. 

This report emphasizes: 

a) How to use DSOP 

b) Block diagrams and listings for better understanding of DSOP 
so that future improvements can be easily implemented 

c) Results obtained on certain test problems 

The following paragraphs discuss the above mentioned items in 
more detail. The 7094 FORTRAN IV program listings are given in Paragraph 6. 

A list of references, with emphasis on those used in the work in this report, 
is given in Paragraph 7. For theoretical background, the reader is referred to 
these references. Many further references on direct search techniques are given 
in the referenced papers. 

2. MAIN PROGRAM 

The primary functions performed in the MAIN program are: 

a) Input-Output 

b) Selection of Type of Search 

c) Pattern Search Logic 

d) Penalty Function Cycling 

A block diagram of the MAIN program is shown in Figure IV-1. Table IV-1 
is a General Symbol Table which explains the most used symbols. Table IV-2 con¬ 
tains the symbols used in VA04A and such key symbols as U, P etc. The under¬ 
lined symbols are vectors, each with the dimension N, where N is the number 
of independent variables. 

The first operation performed by MAIN is initialization and card 
reading. The formats for the input cards can be obtained from the listing in 
Paragraph 6. V> T ~ 
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Figure IV -1 Block Diagram of Main Program for DSOP 






A IV-4 


GENERAL SYMBOL TABLE 
Table IV-1 


FORTRAN 

Symbol 

Programs 

Where Used 

Definition 

BIAS 

MAIN 

Controls minute improvement exit test. 

BVA 

UNIVAR, 

EXPLOR 

Input specifying if U components are to 
be fixed or varied. 

BVB 

UNIVAR 

Records which U components are at 
bounds. 

BVBS 

MAIN 

Storage vector for BVB. 

DELE 

OBJECT 

Intermediate quantity used in PFV. 

IF IN 

MAIN 

Counter used in minute improvement 
exit test. 

IPENF 

MAIN 

Penalty function cycle counter. 

IPENFK 

MAIN 

Number of penalty function cycles to be 
tried. 

K1 

EXPLOR 

STF Parameter. 1/2N < K1 < 1/N. 

K2 

•EXPLOR 

STE Parameter. 

K3 

EXPLOR 

Success STE Parameter K3 >1. 

K4 

EXPLOR . 

Parameter controlling interpolator/ step 

K5 

MAIN 

Parameter controlling Pattern Move 
step size. 

KI1 

UNIVAR 

STE Parameter. 

LCCONV 

SUCCES, 

MAIN 

Specifies state of convergence for last 
iteration. 

LCPF 

VA04A 

Controls retaining W for multiple cycles, 

LC1 

SUCCES, 

MAIN 

Controls exit on failure to reduce P. 
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Table IV-1 / General Symbol Table (continued) 


FORTRAN 

Symbol 

Programs 
Where Used 

Definition 

LC3 

SUCCES, 

EXPLOR, 

UNI VAR 

Controls exit for converged. 

LC4 

MAIN 

Type of basic search. 

LC5 

OBJECT 

Selects an objective function. 

LC6 

MAIN 

Selects UNIVAR or EMR for Pattern 
Move. 

LC7 

MAIN 

Selects plotting for Pattern Move. 

LC8 

' UNIVAR 

Selects New Block Search. 

LC11 . 


Not used. 

LC12 

UNI VAR 

Controls entry to quadratic fit. 

LC 14 

OBJECT 

Controls initialization for orbit trans¬ 
fer O. F. 

LC 15 

OBJECT 

Controls printout of objective function 
title. 

LC17 

OBJECT 

Controls initialization for Test Problem 
#1. 

MINC 

UNIVAR, 

EXPLOR, 

VA04A 

Minimum magnitude for stepsize in the 

U coordinates. (Also see Table 2). 

NPT 

SUCCES 

Counter for points to be plotted. 

OCOUNT 

OBJECT 

Objective function evaluation count. 

OFV 

OBJECT 

Objective function value using penalty 
functions. 

PA, PB 

UNIVAR, 

EXPLOR, 

MAIN, 

NEWB 

Remembered past yaNes ox P. 



A IV-6 


Table IV-1 General Symbol Table (continued) 


FORTRAN 

Symbol 

Programs 
Where Used 

Definition 

PC 

ONEDS 

Remembered past value of P. 

PE 

MAIN 

Remembered value of PA for Pattern 
Move, 

PF 

ONEDS 

Final iteration value for P. 

PFV 

OBJECT 

Penalty function value. 

PT 

EXPLOR, 

UNIVAR, 

NEWB 

Remembered past value of P. 

PFCON 

OBJECT 

Penalty function parameter. 

STE 

EXPLOR 

Vector used in generating n steps ,f 
which change U. 

STF 

EXPLOR, 

NEWB 

Vector used in generating interpolatory 
step. 

STG 

NEWB 

. 5 (UA - UB) calculation for STF. 

SUM 

EXPLOR 

Used in interpolatory step decision. 

UA 

EXPLOR, 

ONEDS, 

MAIN, 

NEWB 

The "advance-point” - the best point 
found to date during the current ex¬ 
ploration. 

UB 

EXPLOR, 

ONEDS, 

MAIN, 

NEWB 

,f Base point” (from which) exploratory 
moves begin. 

UE 

MAIN 

Remembered values of UA in pattern 
search. 

UF 

ONEDS 

Final iteration values for UA. 
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Table IV-1 General Symbol Table (continued) 


FORTRAN 

Programs 


Symbol 

Where Used 

Definition 

UT 

NEWB 

Used for UB storage. 

UMAX 

UNIVAR 

Input upper bounds for U 

UMIN 

UNIVAR 

Input lower bounds for U 


components. 
components. 


XCOUNT 


MAIN 


Counter for. number of UNIVAR or EMR 
iterations. 
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Table IV - 2 

SYMBOL TABLE FOR VA04A 


FORTRAN Math Statement Sequence Numbers 

Symbol Symbol Definition Where Set in VA04A Listing 


A 


Intermediate qty. Ill 

used in D & DD 


* AAA 


Convergence 187, 192, 200, 20b, 209 

Control Parameter 


B 


Intermediate qty. 
used in D & DD 

112 



COUNT 

MAXIT 

No. of iterations 
allowed 

Fixed, Input, 23 7 



D 

d 

L.S. Increment 
Parameter 

43, 60, 73, 108, 

141, 

172 

DA 

a 

Direction Increment 
from DB to FA 
Coordinate 

47, 71, 96, 100, 

105, 

115, 

DACC 


L.S. Parameter for 
Absolute Accuracy- 
Cutoff 

37 



DB 

b 

Midpoint for Quadratic 
Fits 

66, 69, 92, 107, 

117, 

134 

DC 

c 

Direction Increment 
from DB to FC 
Coordinate 

90, 100 



DD 

D 

Change of D at current 
step & 2 £ ^ derivative 
estimate 

48, 142, 168 



DDMAG 


Used in DMAG 
Calculation 

5, 84, 195, 227, 

235 


DDMAX 


Ten Times DMAG 

80, 83, 86, 103 




(1) From Reference 5 
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FORTRAN 

Symbol 

DI 

DL 

DMAG 

DMAX . 

ESC ALE 

FA 

FB 

FC 

FHOLD 

FI 

FKEEP 

FP 

FPREV 


Table IV-2 Symbol Table for VA04A (continued) 


Math. 

Symbol Definition 

Used in Absolute & 
Relative Accuracy 
Tests 

Last Value of D 

Maximum Allowed 
Step Size 

Used in DMAG 
Calculation 

Limits Max. Change 
of Variables 

f P value at an end 

point in 1. s. 

f P value at midpoint 

b . 

in 1. s. 

f P value at an end 

Q 

point in 1. s. 

Value of F after N 
1. s. & before s. 
in new direction 

Intermediate min. 

P setting 

Value of F at previous 
convergence, used in 
complex convergence 
procedure when 
ICON = 2. 

Value of P at beginning 
of iteration 

Value of P at beginning 
of 1.s. 


Statement Sequence Numbers 
Where Set in VAQ4A Listing 

120, 123 

42, 49, 102 
38, 39 

36 

Fixed Input, 5, 6, 18 

46, 70, 95, 99, 104, 114, 131 

65, 68, 91, 106 , 116, 133 

89, 109 

161 

121, 124 

25, 224, 252 

27, 248, 251 

156, 157 
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Table IV-2 Symbol Table for VA04A (continued) 


FORTRAN 
Symb ol 

Math 

Symbol 

Definition 

Statement Sequence Numbers 
Where Set in VA04A Listing 

ICON 


Controls restart for 
accuracy test 

Fixed Input, 214 

IDIRN 


Index for s. Direction 
Components in W. 

34, 146, 183 

ILINE 


Index for s. Direction 

35,. 199, 197 

IND 


A control integer used 
when ICON - 2 

11, 215, 233 

INN 


A control integer used 
when ICON = 2 

12, 217, 267 

IPRINT 


A parameter con¬ 
trolling printing 

Fixed Input, 150 

IS 


Controls Search Mode 

77, 79, 82, 136, 138, 162, 

45 

IS GRAD 


2 for 1st iteration, 

1 thereafter 

23, 235 

ITERC 


Count of iterations 

22, 236 

IT ONE 


New direction control 

26, 126 , 184 

IX P 


Index for W (N + !) 
components 

29, 31, 163, 165, 186, 189, 
199, 203, 253, 255 

J 


Controls integration 
of new direction 

174 

JIL 

m 

Value of ILINE for 
max. SUM 

158, 247, 250, 264 

JJ 

N 2 + N 

Index for new 
direction components 

7 

JJJ 

JJ + N 

Index for W (N + 2) 
components 

8, 243 
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Table IV-2 Symbol Table for VA04A (continued) 


FORTRAN * 
Symbol 

Math 

Symbol 

Definition 

Statement Sequence Numbers 
Where Set in VA04A Listing 

K 



A parameter for DO 
loops 

9, 19, 50, 53, 178, 185, 193 
218, 220, 256 

MINC 

E 


E(I) specifies 
required accuracy 
of argument 

Fixed Input, 16, 205, 206 

N 

n 


Dimensionality of 
space 

Fixed Input 

NFCC 



Count of Function 
Evaluations 

56, 225 

P 

F = 

f(p) 

P(U) - the function 
to be minimized 

55, 140, 201, 225, 241 

SCER 



f(ESCALE) used in 
calculating DACC 

6, 37 

SUM 

A 


Max [fPREV - p] 

ILINE 1— N 

28, 157 

U 

p 


N - component 
argument vector 

52, 144, 204, 222, 261 

w , A2) 

W (o) 



Holds direction 
scaling variables 

18, 148, 196 

W(j) (2) 

£l* • ‘ 

£n 

N stored search 
directions 
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As shown in Figure IV-1, the program input LC4 selects the search 
procedure to be used. Setting LC4 = 5 selects VA04A, LC4 = 4 selects Accel¬ 
erated Convergence, LC4 = 3 selects PMS, LC4 = 2 selects EXPLOR with 
ONEDS, and LC4 = 1 selects EXPLOR alone. Additional options can easily be 
added. All of these search procedures are subroutines except PMS and Accel¬ 
erated Convergence, which are built into MAIN. The operation of PMS will be 
explained in the next section. PMS is not a self-contained optimization algor¬ 
ithm; it requires an "exploratory" subroutine which performs a restricted 
local search about a given point. The two exploratory subroutines presently 
included in DSOP are EMR (Exploratory Move Routine) and UNIVAR. EMR is 
modeled after the exploratory subroutine used by Hook and Jeeves (Ref. !}* 
UNIVAR is similar to part of a published optimization algorithm called BEST 
UNIVAR (Ref. 3 ). These subroutines are described in Paragraph 4. ,f Accel- 
erated Convergence" was an experimental procedure combining a modified 
UNIVAR (produced by setting program input LC8 positive) with a special sub¬ 
routine called NEWB. The Accelerated Convergence method, which was an 
attempt to incorporate part of Powell 1 s method into UNIVAR, produced gen¬ 
erally favorable results in experiments with test functions (see Table IV-3 and 
Paragraph 5) but has been superseded by VA04A. 

PMS and UNIVAR incorporate logic for enforcing upper and lower 
limits on each component of the argument vector. The other search algorithms 
do not include any provisions for constraints. However, regardless of the 
search algorithm selected, MAIN provides means for handling general con¬ 
straints by use of penalty functions. In the penalty-function method, the ob¬ 
jective function is modified by adding extra terms which depend on the con¬ 
straints. This enables a constrained minimization problem to be solved (ap¬ 
proximately) by a sequence of unconstrained minimizations. Each cycle of the 
sequence consists of an unconstrained minimization of the modified objective 
function, followed by an adjustment of the coefficients in the penalty terms. 

In the usual form of the penalty function method, these coefficients are scale 
factors which are increased each.cycle. DSOP,experiments have been made 
with a more sophisticated scheme, in which the penalty terms include bias 
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Table IV-3 


TABLE OF THE NUMBER OF OBJECTIVE FUNCTION 
EVALUATIONS FOR SEVERAL D.S. SCHEMES AND TEST 

FUNCTIONS 








Powell 1 s 
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Test 



Conver- 
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Out 
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MINC 

Function 
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Point of 
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New 

New 

equals 

equals 

(Starting Value) 

Value 

Recording 

Test 

Block 

Block 

. 005 

.05 
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1 

N.A. (1) 

31. 

mm 

111. 

16. 

5. 

( 0 , 1) 


C. 

212. 

1 

II 

152. 

34. 

13. 

CUBE 

2 

N.A. ^ 

183. 

451. 

247. 

45. 

7. 

(0, 1) 


C. 

368. 

544. 

247. 

120. 

!4. 

Colville's 

6 

N.A. (1) 



1039. 

814. 


T.P. #4 

(-3, -1, -3, -1) 


c. 



1235. 

944. 


Colville's 

7 

N.A. 



555. 

3075. 

333. 

T.P. #1 


C. 




(PFCON= 


( 0 , 0 , 0 , 0 , 1) 

1 cycle only 






1000.) 



N.A. - Near Answer 
C. - Cutoff Point 

(1) - Under .5% error in each component {.0005 for 0. minimum) 

(2) - Objective function value close to correct value but coordinates in 

error (also ICON = 1). 
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coefficients that are adjusted each cycle, in a manner designed to give rapid 
convergence to the constrained minimum. 

Figure IV-1 shows the logic which handles penalty function cycling. 
Each cycle, the hopefully improved values of U and F are stored in UB and 
DELE and the flow is routed to block 64 to start the next cycle. OBJECT is 
called again in block 43 because the last computed penalty function values may 
not correspond to the best obtained values of U. The input IPENFK controls 
the number of cycles. This automatic cycling has only been used on Test 
Problem #1 but it could easily be generalized for use with most objective func¬ 
tions with implicit constraints. 

The mathematical form of the penalty function, calculation for Test 
Problem #1 in OBJECT originally was: 


P = P , + PFCON £ f (E. t AE.) 

o.t. . * i i 

i 


where: f(x) = min (o, x) 


E. 

i 


"10 

I 

u=i 


a.. u. 
U J 



AE. = o for 1st cycle 

l 

and the summation on i goes over i = 3, 4, 5, 6, 9 only since it can easily be 
shown that the others need not be considered. At the end of the convergence 
cycle, set E^ = f(E^) using the converged U values and start the next cycle. 
However, recent tests have indicated that only i = 5 and 6 were active for 
VA04A trials on T.P. #1, and OBJECT reflects this* Manual cycling was used 
with the Penalty Function Test function shown in Figure IV-6. 


3. PATTERN MOVE SEARCH (PMS) 

The basic logic of Pattern Move Search is shown in Figure IV-2, and its 
detailed logic shown in Figure IV-1. The program makes use of an exploratory 
subroutine, which operates in two different modes. The exploratory subroutine 
starts at a given base point UB, and makes exploratory moves, according to 
built-in rules, in a search for a smaller value of P. The output of the exploratory 
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subroutine is an "advance point" UA and function value PA, with PA £ PB if the 
search succeeds, and with PA = PB, UA = UB if it fails. In "Mode Plus", the 
subroutine keeps trying, varying step sizes and/or other search parameters, 
until it either succeeds, or concludes that PMS has converged to a solution. 

In "Mode Minus", the subroutine stops after one run through its exploration 
sequence, regardless of success or failure. 

Pattern Move Search proceeds by generating a sequence of accepted 
base points, each of which gives a lower function value than the previous one. 

A "pattern move" generates a tentative new base point by displacing the newest 
accepted base point by a vector equal to (or proportional to) the difference be¬ 
tween the newest and next-newest accepted base points. The tentative new base 
point is improved by use of the exploratory subroutine (in "Mode Minus") and 
the result tested to see if it is an improvement. If so, it becomes the newest 
accepted base point, ana the cycle repeats. If it is not an improvement, i.e., 
if the pattern move (together with its associated local exploration) is a failure, 
the program returns to the newest accepted base point (UE) and executes the 
exploratory subroutine in "Mode Plus". If this fails, the search ends. If it 
succeeds, the output and input of the subroutine become the newest and next- 
newest accepted base points respectively. 

A special merit of Pattern Move Search is that the tentative new 
base point generated by a pattern move is not tested for success until after an 
attempt has been made to improve it by a local exploration. This means that 
Pattern Move Search can follow a curving valley; even if a pattern move misses 
the valley floor and hence gives a high function value, the subsequent explora¬ 
tion sequence may find the valley again, so the over-all move is a success. 
Search programs which always reject points that do not show an immediate 
improvement will sometimes be much less efficient than Pattern Move Search. 
(Sometimes they may be more efficient; persistence is not always a virtue.) 

Most of the logic in. Figure IV-1 can be understood from the pre¬ 
vious discussion, and comparing with Figure IV-2. FORTRAN reference num¬ 
bers to the left of the blocks are used in the following discussion. LCCONV is 
a branching control quantity set by the exploratory subroutine in subroutine 
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SUCCES, and subsequently tested by PMS to decide whether to continue or to 
terminate the search. Another possible termination of the search is by the 
Minute Improvement Test, which counts the number of consecutive iterations 
that have failed to improve the function value by amounts above a given thresh¬ 
old, named BIAS. A third mode of termination occurs if the number of explora¬ 
tory moves, XCOUNT, exceeds a given limit, COUNT. 

Block 17 is initialization* Block 27 saves the newest accepted base 
point (labeling it UE ) and performs a pattern move. The coefficient Kg in Block 
27 is currently set at Kg = 2. By choosing Kg greater than 2, the lengths of a 
sequence of pattern moves can be made to grow more rapidly, but with increased 
probability of failure. Block 22 restores the saved base point after failure of 
the pattern move. The vector BVB , which remembers which components of 
the argument vector are at their bounds (and hence requires different treatment 
in the exploratory subroutine) is saved and restored along with UE . The opera¬ 
tion denoted by BOUND U, in Blocks 55 and 27, is a subroutine that examines 
each component of U to see if it exceeds the prescribed upper or lower bound. 

If so, U(I) is set equal to whichever bound was violated, and BVB is changed 
accordingly. 

Pattern Move Search provides a very general framework into which 
almost any type of optimization algorithm can be fitted as an exploratory sub¬ 
routine, v/ith minimal changes in PMS itself. For example, if it should be de¬ 
cided to modify the optimizer so it can handle general sets of equality and/or 
inequality constraints, this could be done by developing a new exploratory sub¬ 
routine with the desired capabilities. The only changes required in PMS itself 
would be (a) replacing BOUNDU with a more elaborate operation which enforces, 
or approximately enforces, all active constraints, by projecting the search 
point onto constraint surfaces, and (b) saving and restoring information that 
remembers which constraints are active; this requires an extension of the 
BVB vector. 

The PMS program and its associated subroutines were first de¬ 
signed for minimization without constraints, and modified dntii they accorri- 
plished this reliably and moderately efficiently on a set of test problems. 
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From these experiments, UNIVAR appeared to be a better ,r valley-follower" 
than EMR and was therefore selected for further development. PMS and UNIVAR 
were modified to provide for the enforcement of specified upper and lower bounds 
for each component of the argument vector. Also, PMS was modified to incor¬ 
porate an improved stopping rule (the "Minute Improvement Test"). EMR was 
not modified, but EMR and the original version of UNIVAR are still usable for 
problems where argument bounds are absent, or do not affect the solution. 

4. SUBROUTINES 

The ten subroutines and their interrelationships are indicated in Fig¬ 
ure IV-3. The arrows indicate which subroutines are called by each individual 
subroutine. The subroutines are divided into two groups: the primary* search 
techniques, and accessory subroutines which accomplish smaller tasks. The 
subroutine OBJECT evaluates the objective function, and therefore is used by 
all the search programs. Block diagrams of the key subroutines are included 
in Figures IV-4. IV-5 and IV-7. 

4. 1 VA04A (POWELL’S METHOD) 

Powell T s method is what is known as a "conjugate direction" method. 

It proceeds by a sequence of one-dimensional searches. The search directions 
% 

are chosen in a manner which will find the exact minimum of a quadratic func¬ 
tion of n variables with continuous second derivatives, in a finite number of 
steps. For functions which are not quadratic, but can be approximated by- 
quadratic functions in the neighborhood of the minimum, rapid convergence 
may be expected. VA04A has no provisions for handling constraints except by 
penalty functions. * 

A detailed flow diagram of VA04A, as defined by Powell’s unpub¬ 
lished program listing (Ref. 4) is given in Figure IY-4. The principal symbols 
are listed in Table IV-2. To understand the program logic, which is fairly com 
plex,Figure IV-4 and Table IV-2 should be studied together with Powell's paper 
(Ref. 2 ) which explains the principles behind his program. The basic operating 
cycle of VA04A has been concisely described in a review paper by Fletcher 
(Ref. 5) who presents it in pseudo - ALGOL notation as: 
y: = x; 

for i: = 1 step 1 until N do MIN(i); 
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for i: = 1 step 1 until N-l do p^: = p. + 1; 

p N : = y = x; MIN(N). 

where x is an N-component argument vector (equivalent to the U of our notation), 
y is the value x had at the beginning of the cycle, and MIN(i) is a subroutine 
which performs a one-dimensional minimizing search parallel to p^, starting 
from the best point found previously. The vectors p^ are originally chosen to 
be parallel to the coordinate directions. The routine performs a linear search 
parallel to each p^ in sequence. The net resulting change in x defines a new 
p-vector, and all the p-vectors shift down in the list, the oldest (the previous 
p-) being discarded. 

Powell*s actual procedure is somev/hat more complex than this. It 
is evident that if his set of vectors p^ ever become linearly dependent, his 
search process becomes trapped on a hyperplane and will never find the solu¬ 
tion if it is not on this plane. To avoid this, he does not always discard the 
oldest vector; another vector, chosen by a rather complex set of rules, is 
sometimes discarded instead. Also, the newly generated direction is not al¬ 
ways accepted. 

Powell*s program performs N+l one-dimensional searches per 
iteration, each search requiring at least two new function evaluations. It will 
find the exact minimum of a quadratic function in N iterations, requiring a 
total of + N linear searches. For non-quadratic functions, the exact answer 
is generally not found, and more than N iterations will usually be required to 
obtain the accuracy desired. The convergence is ultimately quadratic. That 
is, if the function to be minimized behaves like a quadratic function in the 
neighborhood of the optimum, then if the search point is sufficiently close to 
the solution and at least N cycles have occurred, each subsequent iteration 
will double the number of significant figures. The number of linear searches 
required for minimizing a quadratic function could be halved by a slight change 
in the program logic, but this change would not double the program*s efficiency 
for minimizing non-quadratic functions; at most, it might save half of the first 
+ N searches that would otherwise b oper for me d, but subsequent operation 
would not be improved. 
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The version of VA04A used in DSOP differs from Powell’s original 
program by two small changes. The first of these changes eliminates a diffi¬ 
culty that would otherwise occur in unusual conditions: The zero output of the 
branch on FP-P, at Block 96, is now treated like a negative output (GO to 37). 
This change prevents a division by zero which would otherwise cause the program 
to fail. The second change was introduced to facilitate the efficient use of 
VA04A with penalty-function cycling. It provides an option of bypassing Powell’s 
initialization, so search directions and scaling are preserved from the last use 
of the YA04A routine. 

4.2 UNIYAR 

The subroutine UNIVAR exists in two versions: an original version 
which has no provision for bounds on the argument variables, and a revised 
version that incorporates such bounds. Either version works with PMS, but 
only the revised version, whose flow diagram is shown in Figure IV-5, is included 
here. 

The original UNIVAR proceeds by making a sequence of one-dimen¬ 
sional searches, each parallel to a coordinate direction, until all directions 
have been tried except those for which BVA(I) = 0. Each one-dimensional 
search begins from the best point found by the previous search. The logic for 
a one-dimensional search is as follows: 

Step in the stored direction (and by the stored amount; both are de- 

/' y 

fined by a component of the step vector STE ). If this succeeds, increase step 
size and repeat. If it fails, reverse direction and go the other way. Proceed 
until the minimum along the search direction has been bracketed, i.e., until 
a set of three points has been found with the smallest function-value in the 
middle. This occurs when a success is followed by a failure, or when the first 
step fails and is immediately followed by a second failure when the reverse di¬ 
rection is tried. Fit a quadratic function to the values found at these three 
points, and try the point where this quadratic is minimal. The search result 
is either this new point or the previous middle point* whichever has the least 
function-value. The step size Is reduced t6 cancel the last increase and is 
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further reduced by an extra factor of one-half if the last trial (from the quad¬ 
ratic fit) is not a success. However, step size is not allowed to go below a 
prescribed minimum size. 

This one-dimensional search procedure is known to be inefficient, 
since it requires a minimum of three new function evaluations per direction, 
whereas Powell ! s one-dimensional search procedure (based on remembering 
and updating an estimate of the second derivative) reduces the minimum to two. 
Adopting PowelPs rules for one-dimensional searches might reduce the number 
of function evaluations for each execution of UNIVAR by as much as 33%. The 
average reduction would be less than this, but still significant. This improve¬ 
ment has not been incorporated because other developments have had higher 
priority. 

The revised UNIVAR is the same as the original UNIVAR except 
for extra logic to handle upper and low r er bounds on the independent variables. 
The binary-valued vector BVR keeps track of which variables are at their 
bounds. If BVB(I) = 0, indicating that the Ith coordinate is in the interior of 
its allowed range, the search along the Ith direction is the same as in the orig¬ 
inal UNIVAR, unless a tentative step reaches or exceeds a boundary. In this 
case, a step to the boundary is tried instead. If this step fails, the search 
along the Ith direction will end at an interior point. If it succeeds, BVB(I) is 
set to 1, STE(I) is set to call for a small step (4 times the minimum size) away 
from the boundary, and the routine goes on to the next direction. 

If BVB(I) = 1, indicating that the Ith variable is at a boundary, this 
coordinate is left unaltered unless LCI = +, in which case a step away from the 
boundary is attempted. If this succeeds, BVB(I) is set to 0 and the search pro¬ 
ceeds as from an interior point. If it fails, step size is reduced to half its pre¬ 
vious magnitude, or to the minimum magnitude allowed, whichever is greater. 

The motivation for restricting attempted moves away from a bound¬ 
ary to cases when LCI - + (i.e., to times when the pattern-move sequence is 
being started for the first time, or is being restarted from an old base point 
after failure of a pattern move) is to reduce the number of unsuccessful attempts 
to move away from the boundary that would otherwise occur if the true final 
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solutions call for one or more variables to be at their bounds. Since search 
termination on step size is permitted only when LC1 = +, the program cannot 
exit on step size without having tried moves away from all boundaries. It can 
"exit on test n if there is a prolonged sequence of pattern moves which are all 
successful but give very small improvement, but this is highly unlikely; in gen¬ 
eral, one of the pattern moves will fail, and the next execution of UNIVAR will 
attempt moves away from the boundaries. 

4.3 EXPLOR 

This was the first fT exploratory n program tried with the Pattern 
Move Search. Its principles are described in Refer ence 1. No exercised constraints 
are allowed. Test runs generally indicated EXPLOR was inferior to UNIYAR 
{see Paragraph 5 of this Appendix). The EXPLOR listing in Paragraph 6 can be 
used to follow the discussion below. 

Components of the argument vector U for which BVA(I) = 0 are 
treated as constant constraints; other components of U are varied systematically 
in a search for a local minimum. For each value of I for which BVA(I) •*= 1, the 
routine tentatively changes U(I) in one direction and then (if this fails to reduce 
the function value) in the opposite direction. The direction tried first is deter¬ 
mined by the sign of STE(I), and the magnitude of this quantity determines the 
step size. If either change succeeds, the value of PA (which is the best value 
of P found so far) is updated, and the search proceeds from the altered argu¬ 
ment vector, by considering the next value of I. If both directions of change 
along the Ith coordinate axis are unsuccessful, no change is made in U(I), and 
the routine moves on to the next value of I after saving some auxiliary infor¬ 
mation. 

When all variable components of U(I) have been treated as described 
above, the routine decides whether to accept the result, or try one additional 
"interpolatory” step, whose components are given by a vector STF. The value 
of STF(I) is zero for any value of I for which BVA(I) = 0 or for which either of 
the steps iSTE(I) succeeded. For the remaining cases (i.e., BVA(I) = 1, and 
± STE(I) both unsuccessful), STF(I) is computed by considering the loeatiorr -ofbir. crv 
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the minimum of a quadratic function fitted to the function values found at the 
three points U(I) and U(I) +_ STE(I), with all other components of U fixed. If 
Kq were cho*sen to be 1/2, and no other coordinates were involved, STF(I) 
would give a move to this minimum. A value of Kq less than 1/2 is appropriate 
in multidimensional cases. If P(U) is a quadratic function, and if all coordinate 
steps have failed, it can be shown that the interpolatory step alwa.ys succeeds 
if Kj < 1/N, and is always too short (undershoots the minimum) if Kq < 1/2N. 
Therefore, Kq should be between 1/N and 1/2N. 

The decision to try the interpolatory step, or not try it, is made by 
comparing the total function reduction obtained by steps jb STE(I) with the ex¬ 
pected further reduction obtainable by the interpolatory step. This expected 
reduction is proportional to SUM, if interaction between coordinates is ignored. 
A coefficient Kq is provided which can be adjusted to bias the decision. If 
Kq < 0, the interpolatory step is never tried unless all the steps ± STE(I) failed. 
If Kq is made very large, the interpolatory step will always be tried. 

The step vector STE is updated during the routine, to provide a 
favorable step for next time. If step STE(I) is successful, the new STE(I) is 
larger (K^ > 1). If STE(I) fails but the reverse step succeeds, STE(I) reverses 
sign but does not change magnitude. If both steps +, STE(I) were unsuccessful 
the new STE(I) has its sign determined by comparison of the two unsuccessful 
trials. 

An input option (LC1) is tested in SUCCES to allow the routine to 
exit on failure to reduce the function, or try again with reduced step sizes until 
either success is attained, or the step sizes are reduced to the allowed mini- 
mums given by MINC . 

The performance of EMR could probably be improved significantly 
by adding a one-dimensional search along the direction defined by STE , i.e., 
parallel to the n interpolatory move." However, this experiment has not been 
tried. 

4.4 ONE PS 

.•*. q • ' ...ThjL^ois^.simply • a one-dimensional search program which is used by 

New Block and with a search variation with EXPLOR in the MAIN program. 
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4.5 OBJECT 

This program computes a single number which corresponds to the 

9 

value of the objective function. The inputs are essentially just the independent 
variables U. Seven different test functions (see Figure IV-6) have been used in 
the experimental tests of the optimization algorithms. The selection is con¬ 
trolled by the input integer LC5. The two dimensional version of CUBE was 
taken from Reference 3; Colville T s test problems came from Reference 6 ; the 
orbit transfer problem came from Reference 7« 

Additional functions can be readily added to OBJECT by using new 
L»C7 branches. Also, any available function-evaluation program can be substi¬ 
tuted for OBJECT by appropriately renaming variables and adding COMMON 
statements. 

4.6 NEWB 

NEWS (New Block) uses some of PowelPs basic ideas to accelerate 
convergence. As is shown in Table IV-3, it has been used both with the separate 
accelerated convergence test and in UNIYAR with PMS. 

4. 7 SUCCES 

This program is used by UNIYAR and EXPLOR after every iteration 
to test the condition of the search and set LCCONV accordingly. LCCONV is 
then tested in the PMS in the MAIN program. Also the iteration counter, 
XCOUNT, is updated and information for plotting is stored. See Figure IV-7 for 
the Block Diagram. 

4.8 SETQ 

SETQ performs some simple necessary settings for UA , UB , PA 

and PB. 

4.9 BOUNDU 

This program limits all components of IJ to their bounds, called 
UMAX and UMIN . These quantities must be inputted on cards (or by modifica¬ 
tion through COMMON). BOUNDU is called Exclusively from MAIN, initially 
and in PMS. -Similar operations are also performed in UNIYAR. 
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Table IV-4 

A. R. COLVILLE'S FORMAT FOR 
TEST PROBLEMS # 1 AND # 4 


Test 

Problems & 
Search Method 

MINC 

OBJ, 

Function 

Value 

Prep. 
Time 
Guess 
(Hr s.) 

Execute 

Time 

(Sec.) 

Funct. & 
Constr. 
Calcs. 

Std.< 2 > 

Time 

Ratio 

Test Problem #4 

Powell 

(ICON = 2) 

.os' 1 ’ 

0 . 

.5 

1.3167 

330. 

.0196 

Test Problem #4 
Pattern with 

UNI VAR 
(ICON = 2) 

. 005 

0 . 

.5 

3.667 

1235. 

. 0545 

Test Problem #1 
Powell (1 cycle) 
(ICON = 1) 

.os' 1 ’ 

-32.38 

1 . 

2.450 

333. 

. 0364 


1 MINC = .05 for Powell is equivalent to .005 for Pattern Move since Powell 
tests for . 1 (MINC = E) 


2 7094 Standard Time 


67.2 sec. 
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Table IV-5 

ACCURACY AND TIMING COMPARISONS BETWEEN 
UNIVAR AND EXPLOR WITH PMS 

Test Problem = CUBE 


Starting Conditions 


Running Time (sec.) 
UNIVAR EMR 


Case 1 (0, 0, 0) 
Case 2 (10, 0, 0) 


2.583 

1.683 


2.317 

2.466 


UNIVAR 

EXPLORATORY 

oco 




OCO 




UNT 

P 

U(l) 

U(2) 

UNT 

P 

U(l) 

U(2) 

7 

1. 

0. 

0. 

6 

0.9999 

’ A - 4 

0. 1056. lU 

0. 

1 

1. 

0. 

0. 

18 

0.9996 

0.2082. 10 ~ 3 

0. 14021. 10 ~ ‘ 

31 

0.4727 

0.31250 

0.03052 

30 

0.9905 

0.4779. 10 ~ 2 

0.10216. 10 ~ 7 

51 

0.4172 

0.35408 

0.04439 

54 

0.5849 

0.23575 

0.01031 

79 

0. 1901 

0.56310 

0.17855 

78 

0.2758 

0.47495 

0.10603 

113 

0.07818 

0.72039 

0.37386 

112 

0.06079 

0.75354 

0.42716 

147 

0.02546 

0.84044 

0.59365 

144 

0. 00223 

0.95286 

0.86495 

180 

0.00394 

0.93726 

0.82334 

177 

0.606.^ ® 

1.0075 

1.02294 

220 

0.655. 10-6 

0.99919 

0.99757 

199 

0.570. 10 " 4 

1.0075 

1.02281 

8 

1. 

.596. 10 " 7 

0. 

10 

0. 57554 

0.26401 

0. 

23 

1. 

0. 

0. 

24 

1.4454 

0.26401 

0. 

35 

1. 

0. 

0.254. 10 " 8 

34 

0.6165 

0.26401 

0. 

55 

0.47266 

0.3125 

0.030518 

55 

0.3370 

0.42102 

0.078811 

82 

0.28749 

0.46381 

0.09978 

84 

0.3008 

0.45322 

0.088791 

114 

0.12412 

0.64769 

0.27170 

115 

0.1194 

0.65488 

0.279237 

147 

0.04666 

0.78398 

0.48186 

149 

0. 0255 

0.84058 

0.59484 

182 

0.01184 

0.89119 

0.70781 

183 

0.74.10’ 5 

1.0008 

1.0028 

238 

0.655.10 

0.99919 

0.99757 

204 

0.11.10- 5 

1.0009 

1.0026 


Case 

2 
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4. 10 PLOTOP 

UA, UB , PA and PB from UNIVAR and EMR can be plotted as func¬ 
tions of XCOUNT by setting LC7 + . The necessary storage of quantities is done 
in SUCCES. 

5. COMPARATIVE TABLES 

Four of the test functions shown in Figure IV-6 were used with var¬ 
ious search methods to determine the number of objective functional evaluations 
required since this factor is important for problems with the complicated ob¬ 
jective functions expected to be encountered in the QRGTS. Observations from, 
the results shown in Table IV-3 include: 

1) Powell's method was superior in almost all unconstrained 
cases, as was expected from a review of References 2 and 8, 

2) Accelerated Convergence approached the minimum faster than 
PMS but PMS terminated sooner. 

3) PMS with New Block was generally worse than without it. 

4) Powell's rate of convergence is definitely a function of the 
MINC and ESCALE settings. It appears that MINC in addition 
to ESCALE should be set to rea-sonably high values. 

5) PMS was generally better than VA04A on Colville's Test 
Problem #1 (a problem with constraints) since VA04A failed 
for MINC = 0.5 and ICON = 2 and had an excessive number of 
evaluations for MINC = . 005. 

Table IV-4 is included so Standard Time Ratios of several of th<=> 
runs in Table IV - 3 can be compared with the results in Colville's Nonlinear 
Programming Study summary. All three rank in the top half of his listings 
with .0196 ranking 5th out of 24 methods. The Standard Time Ratio is an index 
developed by Colville (Ref. 6) to compare optimizer efficiencies independently 
of computer speeds. 

Table IV-5 is comparisons between UNIVAR and EMR when used 
with PMS on CUBE. These results, along with another"run not shown, were 
the determining factors in the decision to further develop UNIVAR rather than 
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EMR. For CUBE, EMR has lower OCOUNTS but has a U(2) error of 2% in 
Case 1, and its running time is longer on the average. 



A IV 34 



6. FORTRAN IV LISTINGS 

B 

MATNOP - EFN SOURCE STATEMENT - IFNIS) - 



C BOB COFER 3/?2/67 
C MAIN PROGRA“ FOR DSOP 

DIMENSION UEI1C)»BVBS(20S 

COMMON N.UA(IO),UB(lO),PA»P0 t UllO),P,PF,UFl10),LCII 
C » n VA(10)»K1,K2,K3,K4,COUNT,STE(101,MINC,LC l,LC3 

C ,LCCONV,XCOUNT ,E,LC17,LC9,STF{20>,STG(20) , PFCON 

COMMON /FUNCT/ LC5,LC14.OCOUNT,BI AS , NX , LC15 

COMMON /UN I V/ Kll.UCOUNT,UMlN(20 I.UMAXI20J,BVB(20), ILIMI, ILIM 

COMMON /VAP/ E SCALE* I. PRINT »IC ON fLCPF 

COMMON /PLOT/ USI20,500), NPT,PAS<500),PBS(500) 

COMMON /PE N F/ EX(IC),DELE(10),FUO).OFV.PFV 
REAL MINCn0I,Kl,K2,K3,K4,Kll,NX(10!»K5 
. I FORMAT (2413) 

2 FORMAT (9E 8.51 

3 FORMAT (//(5E20.8)) 

4 FORMAT (IH ///32H0 LIST CF INPUT PARAMETER VALUES) 

5 FORMAT (IK ///17H0 TI ME(SECONDS> =,l£2Q.8) 

6 FORMAT (IH1///43H0 DIRECT SEARCH OPTIMIZATION PROGRAM (OSOP) ) 

7 FORMAT (IH //25H0 LOGICAL CHCICE SETTINGS) 

8 FORMAT (IH /34H0 LC1 (NEGATIVE IS FAILURE EXIT) = ,(1) 

81 FORMAT (42H0 LC4 (TYPE OF BASIC SEARCH - PATTERN-3) *,11» 

82 FURMAT I39H0 LC5 (SELECTS AN OBJECTIVE FUNCTION) * * 11) 

63 FORMAT (45H0 LC6 (SELECTS UNIVAR12) OR EMR WHEN LC4 = 3) =*, III 

84 FORMAT (IH /1I1H0 N UB MINC 

C STE BVA UMIN UMAX ) 

85 FORMAT (7E16.8) 

86 FORMAT (IH //7H Kl -E15.8.7H K2 *E15.8,7H K3 »E15.8,7H K4 
C *E 1 5. 8* 7 H K5 -E15.8/IX.8H COUNT * E15.8.7H BIAS =E15.8,7K Kll 
C«E15.8,9M ESCALE =E15.8,8H PFCCN »E15.8) 

87 FORMAT (IH /26H0 THE INITIAL VALUE OF P «,£15.8) 

88 FORMAT (IH ///41H0 BEGINNING CF DIRECT SEARCH CALCULATIONS) 

89 FORMAT (33H0 SEARCH CONVERGED - BIAS CONTROL) 

91 FORMAT (36H0 SEARCH TERMINATED - COUNT EXCEEDED) 

92 FORMAT I31H0 LC7 (POSITIVE FOR PLOTTING) = ,!1) 

93 FORMAT (32H0 LC8 (POSITIVE FOR NEW BLOCK) *,I1) 

94 FORMAT (33H0 (PRINT (POWELL PRINT CONTROL) = ,I1) 

95 FORMAT I33H0 ICON (POWELL RESTART CONTROL! = ,Il) 

96 FORMAT (46H0 IPENFK INUMBER OF PENALTY FUNCTION CYCLES) « t I2) 

11 RTIM = CLOCKFSOOQ) 

XCOUNT * 0. 

LCCONV * 0 
UCOUNT => 0 
OCOUNT = 0. 

LCll = n 
LC14 * M 
LCl5 * ? 

LC17 * n 
ILIM = r 
NPT * <* 

IPENF = 0 
LCPF * n 

REAO (4,1) N,LC1,LC4,LC5,LC6 »LC7,LC8,!PRINT,I CON,IPENFK 
READ (5,2) Kl,K2,K3,K4,K5 ,COUNT,BiAS *Kl1,ESCALE,PFCCN 
DO 99 I * 1,N 

99 REAO(5,?) NX(I).UB(I),MINC11),STEU),BVA(I),UMIN(I),UMAX{ I) 




26 
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B 

HA I NOP - EFN SOURCE STATEMENT - IFN(S) 


WRITE 

16,6) 


35 

WRITE 

16,4) 


36 

WRITE 

(6,7) 


37 

WRITE 

(6,8) 

LCI 

38 

WRITE 

(6,81) 

LC4 

39 

WRITE 

(6,82) 

LC5 

40 

WRI TE 

(6,83) 

LC6 

41 

WRITE 

(6,92) 

LC7 

42 

WRI TE 

(6,93) 

LC 8 

43 

WRITE 

(6,94) 

IPRINT 

44 

WRI TE 

(6,96) 

ICON 

45 

WRITE 

(6,96) 

IPENFK 

46 

WRITE 

(6,36) 

K1,K2,K3,K4,K5 ,C0UN7,81 AS,K11,ESCALE,PFCON 

47 

> WRITE 

(6,84) 


46 

WRITE 

(6,85) 

(NX f I ) ,UB(I) , MINCH) ,$TEU) , 8VA(I),UHIN( I )«UMAX( I), 


C I * 1,N) 


49 


NC * N «■ 1 
DO 55 T * NCtlC 
55 U(I) * 0, 

C INITIALIZE *n 



IF (LC5 .£0. 3) CALL OBJECT 

00 10 l * 1,N 

71 

10 

um * UBIT) 



CALL BOUNDU 

DO 9 I » 1,N 

81 

9 

ub( ii * um 



CALL OBJECT 

PB * P 

91 


WRITE (6,88) 

** 93 


WRITE (6,87) PB 

94 

50 

DO 51 I * 1, N 


51 

STF(I) * MINCH) 

IF (LC4 # fcO. 3) GO TO 17 

IF (LC4 .EQ. 4) GO TO 52 

IF (LC4 .EQ. 5) GO TO 61 


12 

CALL E*»LGR 

IF (LC4 *EQ* 2) GO TO 15 

PB * PA 

00 13 I * l t N 

114 

13 

UB( I) a UAH) 


14 

XCOUNT a XCOUNT ♦ i. 



IF nxCOUNT .GE. COUNT) .OR. (LC3 .LT. OJ)GO TO 4* 

GO TO l? 


C ONE DIMENSIONAL SEARCH 

15 CALL ONFDS 13* 

PB * PF 

00 16 I * 1tN 

16 UBUI * UFU) 

GO TO 14 

C PATTERN SEARCH 

17 DO 18 T * t,N 

18 UE( I ) * UB (I ) 

LC 11 * 0 
IFIN a o 
LCi * +1 


20 
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B 

MAINOP - EFN SOURCE STATEMENT - IFN(S) - 

IF (LC6 . EO* 2) GO TO 24 

CALL EX&LGR 162 

GO TO 25 

24 CALL UNTVAR 165 

25 IF {LCCONV .EQ. 2) GO TO 43 
IF (LCtONV .EO. 3) GO TO 24 

C EXIT TEST FOR SUCCESSIVE TRIALS WITH MINUTE IMPROVEMENT 

LC1 * -1 

21 IF ((PA-PE+B1AS) .LT. 0.) IFIN * 0 
IF (IFIN .GEc 5) GO TO 41 

IFIN * !FIN ♦ 1 
DO 27 l » t,N 
27 BVBSU) * SVBU) 

DO 19 I ■ IfN 

19 UII) * UE«11 ♦ K5MUAU) - UE(m 

CALL BOUNDU 198 

CALL OBJECT 20C 

00 23 I * t,N 
UECI) * UAU) 

23 U8CI5 - UII) 

PE • PA 
PB * P 

IF (LC6 «EQ« 25 GO TO 31 

CALL EXPLOR 215 

GO TO 3* 

31 CALL UNT VAR 216 

32 IF (XCOUNT *GE• COUNT) GO TO 42 
IF (PA *L7. PE) GO TO 21 

PB * PF 
DO 22 I * i,N 
BVBU) * BVBS(I) 

22 UB( I) * UEU) 

GO TO 20 


41 WRITE (6,89) 236 

GO TO 43 

42 WRITE (6,91) 238 

GO TO 44 

C TEST FOR PENALTY FUNCTION ITERATION 

43 IPENF * IPENF ♦ 1 

WRITE (6,3 ) EX.DELE tF tOFV,PFV 241 

CALL OBJECT 245 

WRITE (6,3 ) EX,DELEfF,OFV,PFV 246 

LCPF * 2 
00 62 T * 1,10 
UBCI) = Utl) 

62 DELE(I) * F (I ) 

44 RTIH2 * CLOCKF{OQQ) 262 

TINE * (RT!M2 - RTIM) • .6 

WRITE (6,5) TIME 263 

IF ((IPENF - 1PENFK) •LT. 0) GO TO 64 

IF {LCT .GT. 0) CALL PLOTOP 268 

GO TO 11 

C EXPERIMENT WTTH ACCELERATED CONVERGENCE 
52 CALL UNI VAR 
DO 53 f * ltN 
UB111 * UA til 


272 
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6 

MAINOP - EFN SOURCE STATEMENT IFN(S) - 

IF (LC6 .EO. 2) GO TO 24 



CALL E*°LOR 

GO TO 25 


162 

24 CALL UNTVAR 

25 IF (LCCONV .EQ. 2) GO TO 43 

IF CLCCONV .EO. 3) GO TO 24 

EXIT TEST FOR SUCCESSIVE TRIALS WITH MINUTE IMPROVEMENT 

LC1 = -1 

21 IF CtPA-PE+BIAS) .LT. 0.) IFJN * 0 


165 

IF CIFIN .GE. 5) GO TO 41 

IFI M = TF IN ♦ 1 

DO 27 T * 1,N 

27 8VS $ ( I) * QVB Cl) 

00 19 l * 1»N 

19 um * DEin ♦ K5i(tAin - uE(m 



CALL BOUNDU 


198 

CALL OBJECT 

DO 23 I * l t N 

UEm * UAUI 

23 UB(I) * um 

PE * PA 

PB * P 

IF CLC6 .EQ. 21 GO TO 31 


20C 

CALL EXPLOR 

GO TO 3? 


215 

31 CALL UNI VAR 

32 IF (XCOUNT .GE. COUNT) GO TO 42 

IF (PA .LT. PE) GO TO 21 

PB * PF 

DO 22 f ■ l f N 

BVBCI) * BVBS(I) 

22 UB< K! * UEU) 

GO TO 20 


218 

41 WRITE (4,89) 

GO TO 43 


236 

42 WRITE (6,91) 

GO TO 44 

TEST FOR PENALTY FUNCTION ITERATION 

43 IPENF * IPENF ♦ 1 


238 

WRITE (6,3 ) EX f DELE»F,CFVfPFV 


241 

CALL OBJECT 


245 

WRITE (4,3 ) EX,DELE,F,CFV,PFV 

LCPF * 2 

00 62 T * 1,10 

uBin * um 

62 DELEU) » F(U 


246 

44 RTIM2 * CLOCKF(QQQ) 

TIME * (RTTH2 - RTIM) • .6 


262 

WRITE (6,5) TIME 

IF ((IPENF - IPENFK) .LT. 0) GO TO 64 


263 

IF (1X7 .GT. 0) CALL PLOTOP 

GO TO 1! 

EXPERIMENT WITH ACCELERATED CONVERGENCE 


268 

52 CALL UNI VAR 


272 

00 53 I * i,N 

U8( I 1 * UA (I ) 





B 

MAIMOP - EFN SOURCE STATEMENT - IFN(S) r 

53 STFCI) * STG(I) 

PB * PA 

IF CtCCONV .EQ.2) GO TO 43 
GO TO 5? 

C POWELL SUBROUTINE 
61 CALL VA04A 
GO TO 4-* 

END 


6 

NEWB1 - EFN SOURCE STATEMENT - IFN(S) 


C MODIFIED POWELL. APPROACH H. ROBBINS. R. COFER 
SUBROUTINE NEfcS 
DIMENSION UT \20) 

COMMON N.UA (10).US (10) .PA .PB tUU 0) « P.PF >UF( 10 ) .LCI 1 
C .BVA(10).Ml,K2.X3.K4.COUNT,STE<10).MINC.LC1.LC3 

C .LCCONV.ACOUNT.E.LC17.LCD *STF{20).STGI20) 

COMMON /UNI V/ KU .UCGUNT«UMIN<20).UMAX(20)•£V0{20>. ILIMI. 1L IM 
REAL K! l.MINCCl 0) tKl .K2 »«3.K4 
PT * PB 
DO It 1*1.N 
UTCD * UB(I) 

II UCI> * UAII) 4 STF(I) 

CALL OBJECT 

IF iP #GE. PA) GO TO 10 
CALL SETO 
CO TO 9 

10 DO 13 1*1.N 

13 UfifI) * U(I) 

PB * P 

9 CALL ONEDS 
PB * PT 
PA * PF 
DO 14 1*1.N 
UBfI) * UT(I) 

14 UACI) * UFCI) 

CALL ONEDS 

PA * PF 
DO 15 I * 1 *>* 

UA C I ) * UF C I ) 

15 STG(I) * .5*<UAU) - U8U>) 

RETURN 

END 
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B 

VA04 - EFN SOURCE STATEMENT - IFN(S) 


C BOB COFER 3/22/67 

C POWELL PROGRAM CONVERTED TO FORTRAN 4 WITH MINOR CHANGES FOR P.F. 
SUBROUTINE VA04A 
DIMENSION W{40 ),Xf10)#E(10) 

equivalence (xm ,um) ,(Ecn #MiNCCin t if,p) 

COMMON N V UA(10) «UB(10) 9 PA 9 P&tU(10) 9 P v PFfUF(10) fLCU 
C fRVA(IO),Ki,K2,K3,K4,COUNT,STE(10),M INC,LCl,LC3 

COMMON /VAO/ ESCALE,IPRINT,ICCNsLCPF 
REAL MINC(10),K1 ,K2 ,K3 ,K4 t Ki 1 , NX < 10 ),K5 
311 FORMAT (5E16.8) 

IF {LC p F .EQ. 2) GO TO 211 
MAX IT * COUNT 
DDNAG^0.1#ESCALE 
SCER*0.O5/ESCALE 
JJ*N»N^N 

NFCC*1 
IND*l 
INN *1 

DO l I * 1 f N 
00 2 J*l t N 
W(K) * 0. 

IF CI-J)4,3,4 

3 WCK)*ABS «E {I ) 1 
W(l5 ®ESCALE 

4 4-1 

2 CONTINUE 
1 CONTINUE 
ITERC-1 
ISGRAD*? 

GO TO 2!2 

211 IND * 2 

212 FKEEP*ABS (F)+ABS (F) 

5 I TONE *1 
FP=F 
SUM*0. 

IXP*JJ 

DO 6 1=1,N 
IXP*IXP*1 

wuxpi»xcn 

6 CONTINUF 
lDIRN*N*i 
ILINE=l 

7 DMAX*Wf TLINE) 

DACC *0MA X* SCER 
DMAG*AMTNi (DOMAG,Q.l*DMAX) 

DMAG*ANA XItDMAG,20.*OACC) 

DOMAX=10 # «DMAG 

GO TO (20,70,71),1 TONE 
70 DL*0# 

0*DMAG 

FPREV=F 

IS»5 

FA*F 


CYCLES 


VA04A006 
YA04A007 
VA04A008 
VA04A009 
VA04A010 
VA04A011 
VA04A012 
VA04AC13 
VAC4A014 

VA04A016 
VA04A017 
VAO4AO10 
VA04A019 
VA04A020 
VA04A021 
VA04A022 
VA04A023 


VA04A025 

VA04A026 

VA04A027 

VA04AC29 

VA04A030 

VA04A031 

VA04A032 

VA04A033 

VA04A035 

VA04A037 


VA04A040 

VA04A041 

VA04A042 

VA04A043 

VA04A044 

VA04A045 

VA04A046 
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OA = DL 

8 DD=D~DL 
DL = D 

58 K*ID1RN 

DO 9 I *1 , N 

X( I ) *X ( I)*DO®W(K) 

K ♦ 1 

9 * CONTINUE 

CALL OBJECT 
NFCC *NFCC +1 

GO TO (10,11,12* 13,14,96) ,IS 

14 IF{F-FA)15 «16 »24 

16 IF (ABS(D) - DMAX) 17,17*18 

17 D~D«D 
GO TO 8 

18 WRITE(6,19) 

19 FORMA T ( *>X ,44H VA04A MAXIMUM CHANGE DOES NOT ALTER FUNCTION) 
GO TO 2" 

15 FB*F 
DB*D 

GO TO 21 

24 F6*FA 
DB*DA 
FA*F 
DA*D 

21 GO TO (83,23),ISGRAD 
23 D*DB*DS-DA 
IS-1 
GO TO 8 

83 0*0*5® fnA+OB-(FA-FB)/(DA-DB)) 

I S*4 

IF((DA-O)®(D-DB))25,8,8 

25 IS*1 

IF(ABS (D-OB)-DOMAX)8,8,26 

26 D*DB*SIOM (DOMA X,08-DA) 

IS*l 

DDMAX*ODMA X+DDMA X 
DONAG*OOHAG*DONAG 
IF(DOMAX-DNAX)8,8,27 

27 DOMAX*DMA X 
GO TO 8 

13 IF(F-FA)28,23,23 

28 FC*FB 
DC*DS 

29 FB*F 
08*0 
GO TO 

12 IFCF-F8)2Sf28»31 

31 FA*F 
OA^O 

GO TO ^ 

11 IF(F-FR)32,10,10 

32 FA*F8 
DA*08 
GO TO 

71 Ot*U 


VA04A047 
VA04A048 
VA04A049 
VA04A050 
VA04A05i 
VA04A052 
VA04A053 
VA04A054 
VA04A055 
VA04A056 
VA04A057 
VA04A058 

VA04A060 
VA04A061 
VAG4A062 
VA04A063 
VA04A064 
VA04A065 
VA04A066 
VA04A067 
VA04AQ68 
VA04A069 
VA04A070 
VA04A071 
yA04A072 
VA04A073 
VA04A074 
VA04A075 
VA04A076 
VA04A077 

VA04A079 
VA04A080 
VA04A081 
VA04A082 
VA04A083 
VA04A084 
VA04A085 
VA04A086 
VA04AG87 
VA04A083 
VA04A089 
VA04A090 
VA04A091 
VA04A092 
VA04A093 


VA04A096 

VA04A097 


VA04A100 
VA04A101 
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ODNAX*5. 

FA = FP 
DA=- 1. 

FB*FHOLP 

06*0. 

0*1 • 

10 FC=F 
OC = D 

30 A s (D8~nt 5 *(F A-FC) 

B = CDC-DA)MFB-FC) 

IF((A+B) *(DA-DC)533,33,34 

33 FA*f B 

FB = FC 
D9*DC 
GO TO 76 

34 0=0* 5* ( A* (DB + DO+B* CDA+DC ) i / ( A+Bl 
Dl*D8 

F I*FB 

IFCFB-fC)44,44 f 43 

43 DI*DC 
F I«FC 

44 GO TO (86,86,855, l TONE 

85 I TONE*? 

GO TO 4*5 

86 IF CABS (D~DI)~DACC) 41,41,93 

93 IF CABS |£-DI)-0.03«ABS (0)) 41,41,45 

45 IF ((DA-DC)*COC-D55 47,46,46 

46 FA»FB 
OA*DB 
FB*FC 
DB*DC 

GO TO 75 

47 IS* 2 

IF I(08-0)»(D-DC) ) 48,8,8 

48 I S*3 

GO TO 8 
41 F*F I 

0*0I-OL 

00*SQRT ( <DC-0B5*(0C-DA)*<DA-DB>/tA*B> ) 

00 49 I=1 ,N 

XC I $*X( T )+0»U(IDIRN) 

KUDIRN5*DD*VH!DIRN) 

IDIRN* T0 IRN + 1 

49 CON TINUF 

W(ILINE)*W(ILINE)/DO 
ILiN£*TLINE*l 
IF CIPRINT-1)51,50,51 

50 WRITE{F.52)ITERC,NFCC,F,(XC!),I*i,N5 

52 FORMAT !/1X,9HITERATION,I 5,i15,16H FUNCTION VALUES, 
i 10X,3HF *,E2i.t4/(5E24.14)) 

GO TOC 51,53 5,1 PRI NT 

51 GO TO (55,38),IT0N£ 

55 If (FPREV-F-SUK) 94,95,95 
95 $UM*FFREV-F 
JIL*ItrNE 


VA04A103 

VA04A104 

VA04A105 

VA04A106 

VA04A107 

VA04A108 

VA04A110 

VA04A112 

VA04A113 

YA04A115 

VA04A116 

VA04A117 

VA04A118 

VA04A120 
VAC4A121 
YA04A122 

VA04A124 


VA04A127 


VA04A132 

VA04A133 

VA04A134 

VA04A135 

VA04A137 

VAG4A139 

VA04A141 
VAQ4A142 
VA04A143 
VA04A144 
VA04A145 
VA04A146 

VA04A148 
VA04A149 
VA04A150 


VA04A134 


VA04A158 
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94 

IF UDISN-JJ) 7,7,84 


84 

GO TO (92,72),!N0 


92 

FHOLD^F 



I S=6 

VA04A162 


IXP*JJ 

VA04A163 


DO 59 1*1,N 

VA04A164 


IX?*£XP«-1 

VA04A165 


«(IXP)*Xm-WCIXP) 

VA04A166 

59 

CONTINUF 



DD* 1 • 

VA04AI68 


GO TO 58 

YA04A169 

96 

GO TO (112,87),INO 


112 

IF (FP-F) 37*37,91 


91 

D*2.»(F*+F-2.*FHQLD)/(FP-F)*«? 



IF (D*(FP-FHOLD-$U^) **2-SUM) 87,37,37 

VA04A173 

87 




IF (J-JJ) 60,60,61 

VA04A175 

6° 

DO 62 



K*I-N 

VA04A177 


HCK)*H(n 

VA04A176 

62 

CONTINUF 



DO 97 T*J!L,N 

VA04A180 


H(I-1)*W(I) 

VA04A181 

97 

CONTINUF 


61 

!DIRN*TOIRN~N 



ITGNE*3 

VA04A184 


K*I DIRN 

VA04A185 


IXP*JJ 

VA04A186 


AAA*0* 

VA04A187 


DO 65 I *1 ,N 

VA04A138 


IXP-IXP + 1 

VA04A189 


NtK)*W(TXP) 

VA04A190 


IF (AAA-A85 (W(K)/E(I))) 66,67,67 

VA04A191 

66 

AAA*AB S (W(K)/EU)) 


67 

X*M1 


65 

CONTINUF 



DDWAGM• 

VA04A195 


W(N)*ESCALE/AAA 

VA04A196 


ILINE*N 

VA04A197 


GO TO 7 

VA04A198 

37 

IXP*JJ 



AAA«0. 

VA04A200 


F*FHOLH 

VA04A201 


DO 99 I*1, N 

YA04A202 


IXP*IXP+1 

VA04A203 


X(I1*X(f)-W(IXP) 

VA04A204 


IF (AAA»ASS CE(I!)-ABS (W(IXP))) 98,99,99 

VA04A205 

98 

AAA *A8S CWUXP)/ E(I)) 


99 

CONTINUF 



GO TO 7* 

VA04A208 

38 

AAA *AAA* € 1•♦D1) 



GO TO (72,106),I NO 

VA04A210 

72 

IF (1PRTNT-2) 53,50,50 


53 

GO TO ( 109,88),I NO 


109 

IF CAAA-0.1) 89,89,76 


89 

GO TO (20,U6) , I CON 
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116 INO *2 

GO TO (100,101) »I NN 

100 INN *2 
K-J J J 

00 102 1*1,N 
♦ 1 

WCK )*X( I | 

xm*xm + io.» eu) 

102 CONTINUE 
FKEEP*F 
CALL OBJECT 
NFCC*NFCC«-1 

0DMAG*O, 

GO TO 108 

76 IF CF-FP) 35,78,78 
78 WRITE(6» 80) 

80 FORMAT (5X,37HVA04A ACCURACY LIMITED BY ERRORS IN F) 
GO TO ?o 

88 INO *1 

35 DDMAG*0.4*$QRT (FP-F) 

ISGRAO*! 

108 ITERC*!TERC*1 

IF CITERC-MAXIT) 5,5,81 

81 WRITE (6,821 MAXIT 

82 FORMAT!I 5,3OH ITERATIONS COMPLETED BY VA04A) 

IF (F-FKEEP) 20,20,110 

110 F*FKEEP 

00 ill !**,N 

JJJ-JJJ+i 

xm*wc jjj) 

111 CONTINUE 
GO TO 20 

101 J IL *1 
FP-FKEF 0 

IF (F-FKEEP) 105,78,104 

104 J !L*2 
FP = F 
F*FKEEP 

105 IXP*JJ 

00 113 1*1,N 

IXP*IX*»+1 

K*IXP+N 

GO TO (114,115) ,JIL 
114 WtIXP)*VI(K) 

GO TO 11 3 

us wuxp)*xm 
xu )=wm 

113 CONTINUE 
J IL *2 
GO TO 92 

106 IF (AAA-O#1) 20,20,107 
20 RETURN 

107 lNN*i 

GO TO 35 
END 


VA04A216 

VA04A218 
VA04&219 
VAC4A220 
VA04A221 
VA04A222 

VA04A224 

VA04A225 

VA04A226 

VA04A227 

VA04A226 


VA04A232 


VA04A235 

VA04A237 


VA04A240 

VA04A242 

VA04A243 

VA04A244 

VA04A246 

VA04A248 
VA04A249 

VA04A251 

VA04A252 

VA04A254 

VA04A255 

VA04A256 

VA04A257 

VA04A259 

VA04A261 

VA04A263 

VAC4A264 


VA04A268 
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C H. ROBBINS. R. COFER 3/22/67 
SUBROUTINE UNIVAR 

C BEST UNIVAR DIRECT SEARCH PROGRAM - tfITH PROVISIONS FOR M IN—MAX CONSTRAINTS 
COMMON N.UA < 10) .UB C 10) .PA.PS.UUO) .P.PF.UFC10 ).LC1 l 
C «BVA{ X 0)* K1 •K2.K3*K4.COUNT.STE (10). MINC.LC1.LC3 

C .LCCONV.XCOUNT.E.LC17.LC8.STF(20).STG{20) 

COMMON /UNI V/ K 1 1 # UCOUNT ,UMI Nf 20) ,'JMAX (20) .BVB( 20) • IL IMI. IL IN 
REAL Kl J .MINCU 0) .K1 *K2 .K3.K4 
4 FORMAT C//(5E18.8>) 

PA = PB 
UCOUNT = 0 

IF (LC I .GE. 0) GO TO 9 
DO 8 I * l.N 

8 STF(I) = STG(I) 

9 LC3 * -2 

DO II I * l.N 

is uin - uocn 

C START INNER LOOP 

DO 29 I * l.N 
IF (BVA(I) .EG. 

IF (DVDCl) #EQ« 

IF (LCI .LT. 0) 

UT * U<11 

UC I) * U(I) ♦ STE (I ) 

CALL OBJECT 
IF (P .LT. PA) GO TO 


0.)GO TC 29 
0.) GO TO 88 
GO TO 29 


86 


GO TO 99 

86 DVB(l) - 0 

87 PT * pa 
PA * P 

DEL2 * STE(I) 

GO TO 90 

88 LC 12 « 2 
GO TO 89 

90 STE(I) * KI!♦STE <I) 

LCI2 * -2 

89 U? * UCI) 

UCI) * UCI) ♦ STE(I) 

DEL * UMINCI) - UCI) 

IF (DEL .GE. 0.)GO TO 94 
DEL = UMAX!I) - UCI) 

IF (DEL .LE. 0.)GO TO 93 
CALL OBJECT 

IFCCP - PA) .LT. 0.) GO TO 87 

92 DELI » UCI! - UT 
UCII * UT 

IF (LC12 .LT. 0) GO TO 98 
PT * P 

STECI) * —STE CI) 

DEL2 * —DELl 
GO TO 90 

93 UCI) * UMAX(I) 

CT * - MINC(I) 

GO TO 95 

94 UCI) * UMIN(I) 
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DT * MINC(I) 

95 CALL OBJECT 

IF UP-PA) •GE. 0*) GO TO 92 
PA * P 
BVBCI) = 1 
STE|I) * DT*4. 

GO TO 29 

98 PMPA = P - PA 
PTMPA s PT - PA 

E M • 5*UD£L2**2)*PMPA - (DELI *$2) ^PTMP A) )/( DEL24PMP A 40£L1*PTMPA> 

UU) a UU) - E 

CALL OBJECT 

STECI ) * STE <I)/ Ki1 

IF UP-PA) *LT. 0.) GO TO 71 

99 STE<I) = STE(I)/2. 

U( I) « UT 

GO TO 72 

71 PA * P 

72 DEL * ABS(STEU))- MINCU) 

IF (DEL eGT. 0.) LC3 - 2 

IF < DEL »LT* 0.) STECI) *s SIGNCMINCU)•STECI)) 

29 COHTINUE 
C E^D OF INNER LOOP 
DO 31 I * l'N 

31 UACl) * UCl) 

IF CLC8 «GT• 0) GO TO 5i 

32 LCCONV * 0 
CALL SUCCES 

IF f LCCONV #EQ« ♦) GO TC 9 
RETURN 

C PROVISION FOR NEW BLOCK - ACCELERATED CONVERGENCE 
51 CALL NEWS 
GO TO 32 
END 
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C EXPLORATORY MOVE ROUTINE - H* ROBBINS, R. COFER 
SUBROUTINE EXPLOR 
DIMENSION STF(IO) 

COMMON N,UAC 10) fUBUO) *PA,P9*UClO) ,P*PF•UF(l0)•LC1l 
C #B VA{ 1 0 ) « K1 ,K2,K3, K4, COUNT,STEC10) ,MINC,LCl,LC3 

C ,LCCQNV, XCOUNT 

REAL MINCC10),K1,K2,K3,K4 
4 FORMAT C//(5E18«8>) 

11 LC3 ft-l 
SUM * 0* 

PA s PB 

DO 12 I * l.N 
UtI) * UBC!) 

12 S7FCS) » 0. 

C START OF FIGURE I 

DO 29 I * I,N 

IF (BVA(I) «EG« 0.) GC TO 29 
UCI) = UU) 4 STE(I) 

CALL OBJECT 

IF CP - PA) 13,14*14 

13 STECl) * K34STECI) 

15 PA * P 
LC3 * 41 
GC TO 29 

14 PT * P 

UCI) * UCI) - 2• *STEII) 

CALL OBJECT 

IF (P - PA) 16,18,16 

16 STE CI) * —STE ( I ) 

GO TO 15 

18 UCI) = UCI) 4 STECI) 

PMPT * P - PT 

PPPT ft P 4 PT 

PPPTPA * PPPT - 2 • * PA 

STF(I) * Kl* PMPT * STE(I)/ PPPTPA 

OO = SIGN<I*,PMPT) 

IF C PMP T ,EG* 0.) QQ = U 

STECI) ft K24STE Cl) *GG 

SUM * SUM 4 CPMPT*42/PPPTPA) 

OELTA « ABSC STE C I ) ) - MINCCI) 

IF COELTA) 21*21,19 

19 LC3 ft 41 
GO TO 29 

21 STECI) * MINCCI)♦SIGN(1.,STECl)) 

29 CONTINUE 
C END OF INNER LOOP 

¥PITE C 6*4) U,STE,STF,P,PA,PT,SUM 
C 6N0 OF FIGURE 1 

DELTA a PA - PB 4 K4*SUM 
IF (DELTA) 31,33,33 

31 00 32 I * 1,N 

32 UA< I) * UU) 

LCCONV ft 1 
GO TO 41 

33 00 34 I = l.N 




B 

EXPL - EFN SOURCE STATEMENT , - IFNCS) - 

34 U<I) = UCI) + STF(I) 

CALL OBJECT 

IF <P - PA) 35*37*37 

35 PA = P 

DO 36 l = UN 

36 UAC I) = U(I ) 

LCCONV = 1 
GO TO 41 

37 DO 38 I 3 1 * N 

38 UA(I) = U ( I ) — STF{I) 

LCCONV = 0 

41 CALL SUCCES 

IF (LCCONV «NE• 4) RETURN 

43 DO 44 I * 1*N 

44 STE(I) s •54 STE (I) 

GO TO I 1 

END 

C BOB COPER 3/22/67 

C WRITES FINAL CONDITIONS AND STORES PLOTTING INFORMATION 
SUBROUTINE SUCCES 

COMMON N.UA <10) ,UOI10) *PA*PB.U(10).P*PF*UF(10)•LC11 
C *BVA(103*K1,K2.K3•K4 *COUNT .STE {10).M1NC.LC1.LC3 

C ttCCONVyXCOUNT»E 

COMMON /FUNC T/ LC5.LC14.OCOUNT*81 AS•NX . LC l 5 
COMMON /PLOT/ US(20*S00)*NPT*PAS(500)*P03(500) 

REAL K11 • MINC(10) * K1 *K2 *K3* K4 * N X(10) 

1 FORMAT (15H0 SEARCH FAILED) 

2 FORMAT (16H0 SEARCH SUCCESS) 

3 FORMAT {4OH0 SEARCH CONVERGED • PA=PS AND LC3 MINUS) 

4 FORMAT (//(5E18.8)) 

5 FORMAT( 1H ///48H0 END OF ONE COMPLETE PASS THROUGH UNIVAR OR £MR) 

6 FORMAT <7H P =E15.8.7H PA *E15 0 8.7H PB = E15.8.9H XCO 

CUNT =EI5.8.9H OCOUNT =£15.8} 

7 FORMAT C58H0 N U UB 

C STE/(4E16.8)) 

WRITE (6*5) 

IF t PA .LT. PB) GO TO 41 
IF (LC3 .LT. 0) GO TO 39 
LCCONV * 4 

IF (LC1 .GE. 0) GO TO 42 
WRITE (6.1) 

LCCONV * 3 
GO TO 42 
39 WRITE <6.3) 

LCCONV - 2 
GO TO 42 
4t WRITE (6*2) 

LCCONV * l 

42 XCOONT s XCOUNT ♦ 1. 

C STORE FOR PLOTTING 

NPT a NPT ♦ 1 

00 43 I * l.N 

43 LSI l.NPT) 3s UU ) 

PA SINPT) * PA 
PBS(NPT) a PB 

WRITE < 6 * 6)P * PA.PB»XCOUNT.OCOUNT 

WRITE (6*7) (NX <I ) •U(I) * UB <I) * STE(I}* l * l.N) 

R6TU3M 

EWO 
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C ONE DIMENSIONAL SEARCH PROGRAM H# RG38INS• R« CCFER 
SUBROUTINE OWEDS 
I FORMAT ( //(5E20•8)) 

COMMON N*UAeiO)*UB{10>•PA,PB 8 U(10),P*PF t UF<lO) 
REAL K 

10 DO It I * UN 

11 UI I) * UA Cl) ♦ 2«*(UA(I} - UB ( i } ) 

CALL OBJECT 

IF (P #GE« PA) GO TO 9 
CALL SE TO 
CO TO 10 
9 PC * P 

DO 12 I * 1«N 

12 U(t) » UA(I) ♦ {UA { I ) - US(II) 

CALL OBJECT 

IF <P •GE. PA) PC * p 
IF iP .LT* PA) CALL SETQ 
K * • 5$ CPB—PC ) / (PB4-PC—2 • ♦PA ) 

DO 13 I « 1*N 

13 UCI) * UACI) + K*<UA(I) - UB (1) ) 

CALL OBJECT 

IF CP - PAl 16*14«14 
|4 PF * pa 

DO 15 I * 1*N 

15 UFI I) * UA(I) 

GO TO 16 

16 PF * P 

DO 17 I « !*N 

17 UFII) * UII) 

16 WRITE Ctcl) UA•UB« U *UF * PA *PB »P»PF «K 
RETURN 
END 


6 

SETQ1 - EFN SOURCE STATEMENT - IFNlS) 


16/31/44 


C BOB COFEft 3/22^67 

C PERFORMS REPETITIVE SETTINGS OF A AND B QUANTITIES 
SUBROUTINE SETQ 

COMMON W * UA HO) * U3 110) *PA*PB*UOO) *P»PF»UF(10) 

PB * PA 

PA « P 

DO IS I * 1»N 
UB!n * UAII) 

UA I !' * 111) 

ftETLRN 

END 


12 
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C BOB COFER 3/22/61 

C LC5 SELECTS ONE OF EIGHT TEST FUNCTIONS 
SUBROUTINE OBJECT 

DIMENSION EUK5) ,EU2(5) ,R2MRU5) »R2(5) t Rl (5) f VU5 ) ,V2t 5) f Vt 5) 

C ,RICR2(5) 

DIMENSION Wi(5 ) v W 2(5)«E1(5) 9 E2(5)«Ul(5)«U2(5) 

DIMENSION AA(10,5>,BBtiO),CC(5,5),00(5) 9 EE(5) 

COMMON NtUA(lO) V UB(10) tPAtPBtUdOS tPfPFtUFClO) V LC11 
C tBVAdO) tKl ,K2 f K3 , K4,COUNT t ST E (10) , MINCt LCl, LC3 

C f LCCONVt XCCUNT.E »LCl 7 ,LC8.STH20),$TG(20),PFCON 

COMMON /FUNCT/ LC5 # LC14,OCOUNT,31 AS,NX,LCl5 
COMMON /PENF/ EX ( 10) ,DELE ( 103 »F dO) »CFV # PFV 
REAL MU*11(5),12(5),IN,MF,HINCCIOJ 9 NX(10 I 

1 FORMAT I9ER.5) 

2 FORMAT C1H MlHO THE OBJECTIVE FUNCTION * ORBIT TRANSFER) 

3 FORMAT CIH /85H0 THE OBJECTIVE FUNCTION * CUBE * 100 • *( UC 2 5~UU )* 

C*3)*«2 * ll«-Ud))»«2 ♦ 10* »UC3> J ■ 

4 FORMAT (1H /63H0 THE OBJECTIVE FUNCTION »U«-U(l ))**2 ♦ ABS(U(2>- 
CU( U*»3)«IGC.) 

? FORMAT (1H /49H0 THE OBJECTIVE FUNCTION *F(’J) ♦A«G(U) *B»G(U)«*2) 

8 FORMAT flH /48H0 THE OBJECTIVE FUNCTION » Ull)«2 ♦ 50.*U<2)**25 

9 FORMAT U/i 5El8.8i) 

99 FORMAT (5E16.85 

C PRINT OBJECTIVE FUNCTION NAME 
IF ILC15 .EQ. 0) GO TO 10 
LCl 5 * * 

IF (LC5 «EQ« 2) WRITE (6,3) 

IF (LC5 .ECU 4) WRITE (6,4) 

IF (LC5 .BO. 5) WRITE (6,7) 

IF (LC5 .EQ* 1) WRITE (6,8) 

10 OCOUNT * OCOUNT «■ 1. 

IF (LC5 .EQ. 2) GO TO 6 

IF (LC5 .EC. 3) GO TO 21 

IF (LC5 .EO. 4) GO TO 42 

IF (LC5 .EQ. 5) GO TO 51 

IF (LC5 .EQ. 6) GO TO 61 

IF (LC5 .EQ. 7) GO TO 71 
P * U(1)**2 ♦ 50.*U(2)**2 
GO TO 5 

6 P * 10n.#(U(2) - U (i) **3) *«2 ♦ d. U(l))*»2 * 10.»U(3)»*2 

5 IF CIF .GE. PA) .OR. (LC11 .EQ. 0)) GO TO 12 
CALL SETQ 

GO TO n 

12 LC11 * 0 

13 RETURN 

21 IF C LC14 * EO. 25 GC TO 30 
C INITIALIZATION FOR ORBIT TRANSFER OBJECTIVE FUNCTION 
WRITE (6 f 2) 

LC 1 4 * 7 

READ (5,U WWltEEl,Pl ,WW2,E£2,P2,!N 
WRITE (6 f 9) WWl,EEl,Pl f WW2 t ee2,P2,IN 
MU * 14*076539E15 
MF * 52*<U 
FM * 1./5280. 

OR » .0174532925 
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RO * 57.2957795 

WWi = 0R*WW1 

ViW2 = HR*WW2 

Pi * MF*Pi 

P2 - MF*P2 

PAR 3 = SQR T(MU/P1) 

PAR4 * ^QRTlMU/P2) 

IN = DR * IN 

Cl = comin) 

SI * STNUN) 

SWWi = SINIWWi) 

00 41 T = 1,2 
STE(I) = DR* STE(I ) 
U8<I) * DR*UB(!) 

41 MINCH) * MINCU)*DR 
STE(3 5 = MF* S TC(35 
U8<3) * MF*U8(3) 

HIMC C 3} * MF* MlNC{3 5 
0. 

- S! 

Cl 

0 . 

0 . 

w 

EE1*C0SCWWU 
SWVi«Cl*EEl 
SWWi* SI*EEi 
EE?*CQS(WVj2) 
EE2* SI N ( VIW2 5 
0. 


Wl( 1) 

W1C 2) 

Wll 3) 

W2( i) 

W2( 2) 

W2( 3) 

e ic n 

EH 2) 

Ell 3! 

E2( l\ 

E2125 
E2( 3) 

GO TO 35 

C START OF REPETITIVE SECTION 
30 SUi * STNCUim 

uim = cosium j 


UK 2) 
Ulf 3) 
U2( 1) 
U2( 2) 
U2( 3) 


31 


SU1*CI 
SUI*SI 
COS(U(2 ) ) 

SINlU(2)) 

0. 

VAR 1 = ♦ EEUCOSlUm - UUD) 

VAR2 = P2/(1. ♦ EE2*C0SIU(2) - WW2) ) 

VAR 3 = SQR T(HU*U(3) > 

VAR4 * SQRT { MU/U(3)) 

00 31 T * 1*3 
Rim * VAR 1*U1 II) 

R21 I# * V.^2»U2(!) 

EUim * E! (!) ♦ U1 (I) 

EU2U5 » E2U) ♦ U21H 
R2MR iff) * R2U) - Rl(I) 

VI t i) * PARS* f Wl C?)*EU1 (3) 

PAR3*{*U tl)»EUi C 3 ) 

PAR 3* IWI 11) *EU1 12) ~ Wl C 2 ) *£U 1(1)1 
PAR 4* (W2 (2 ) *EU2 (3 ) - W213)*EU2(2)) 
V2( 2) =-PAR4*CW2m*6U2 (3) - H2(3)*EU2(U) 

V 2{3) * P4R4*lW2(l)*eU2(2) - H2(25 *EU2(1)) 

R ICR ? 11)=R1(2)*R2(3) - Ri(3)*R2(2) 


VH2) 
VH 3) 
V2CU 


- Wi(3)*EUi(21) 

- wim*EUKin 


IF N C S ) - 
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a 

08JT - EFN SGURCE STATEMENT - IFN(S) - 

R1CR2C2)*-R1C1)*R2(3) ♦ Rl(3)*R2(l) 

RlCR2m * Rl(l)*R2(2) - R1 <21 *R2 C11 

R1CR2H « $QRT(R1CR2(1)**2 4 R1CR2(2)*«2 ♦ RICR2(3)**2) 

DO 32 I * 1,3 

32 V(l) * VAR 3*R2 MR 1(I)/RICR2M 

DT *ARCns<Ul(l)*U2(l) 4 Ui(2)*U2{2) ♦ U1(3)*U2(3)) 

2 * VAR A* TAN ( DT/2. ) 

* DO 33 I * 1,3 

Iim « V(l) 4 Z*U1 (I ) - VI (!) 

33 12(1) * V2U) ~ (V(I) - Z«U2U)) 

P * SQRT(nm*#2 4 11(2)**2 4 11(3) **2 ) 4 SORT ( 12 ( 1 ) * *2 4 I2(2)*» 
C 2 * I2(3)**2) 

35 RETURN 
C TEST FUNCTION 

42 P * (l. - U(1l)**2 * ABS (U (2 ) - U(l)*«3)* 100. 

tCll ■ 0 
RETURN 

C RENAL TV FUNCTION 

5i F * U(ll ♦ U(2) 

A * .795968 
8 * 50. 

G r U(l)#*2 4 U{2 ) »*2 -1. 

P a F ♦ A*0 4 B*G**2 
RETURN 

61 P » IQ0.#(U(2)-U(l)**2)«*2 4 (l..U(i))t#2 4 90.*(U(45-U(3)*®2) 

C 4 Cl.-UC 3))4«2 4 10.i*((U(2)-1»)4*2 4 (U(41-1.)*#2) 4 19.8*(U(2)~ 
CU)MUC'*)~1.) 

RETURN 

C TEST PROBLEM NO. 1 


71 

IF (LC 1 7 . 
LC17 * 2 

EQ. 2) GO TO 74 


DO 72 I * 

1,10 

72 

READ (5,1) 

(AA(I,J) ,J*1,5) ,88(1) 


DO 73 I « 

1 # 5 

73 

READ (5.1) 

ICC(I,J),J = 1,5) 


READ (5,1) 

DD 


READ (5,1) 

EE 


READ (5,1)DELE(3),DELc(4).DELE(5),DELE(6),DELE(9) 
WRITE (6,99) AA.8B.CC .DD.EE 
C COMPUTE OBJECTIVE FUNCTION VALUE 
74 SCYY * n. 

00 77 J * 1,5 
SCY « 9. 

DO 76 I « 1,5 

76 SCY « CC(ItJ)»um 4 SCY 

77 SCYY « SCY*U(J) 4 SCYY 
EPDY * 5 

DO 78 J « 1,5 

78 EPDY * (HJ)MEE(J) 4 OD ( J ) *U ( J) *#2 ) 4 EPOY 
OFV * SCYY 4 EPDY 

C PENALTY FUNCTION CALCULATION 
DO Si I * 5,6 
ext!) * o. 

DO 79 J * 1,5 

79 EXU) * AA f I , J) «U ( J) 4EX( I) 

81 Fin - FXU) - R8(l) 4 DELE ( I ) 
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B 

08 J T - EFN SOURCE STATEMENT - IFN(S) - 

F(3) * 

F(4) * 0. 

F(9) = 

PFV »PFCQN*(F(3)**2 ♦ F(4)**2 «• F(5)«*2 ♦ F(6)**2 + F(9)#*2) 

p = OFV ♦ PFV 

RETURN 

END 


B 

BOUND - EFN SOURCE STATEMENT - IFNCS) 




C BOB COPER 3/22/67 

C LIMITS EACH COMPONENT OF U TO ITS EXPLICIT CONSTRAINTS 
SUBROUTINE BGUNDU 

COMMON N * UA f 1 0 ) • UB ( i 0 ) * PA • PB * U ( 1 0 ) ,P.?F,UFUO)tLCU 
C *B VA ( 10)« K1* K2 * K3 *K4*COUNT cSTEIIO)#MINC *LCI * LC3 

COMMON /UNI V/ K 11 * UC0UN7 § UNI N (20 ) *UMAX(20) • BVB ( 20 ) * 1L IM I * II. IM 
DO 11 I * t*N 
II BVBCX) * 1. 

DO 4® I * t*N 

DEL * UM IN {I ) — UCI) 

IF CDEL ) 43*46*46 

43 DEL * UMAX<I) - U(I) 

IF <DEL > 45*45*44 

44 BVBtl) *0• 

CO TO 48 

45 U< I) * UMAXU) 

STECI) * —AB S(S TE(I )) 

CO TO 46 

46 UCI) * UMIN(I) 

STE i I> * AB S(STE Cl)) 

48 CONTINUE 
RETURN 
END 
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B 

PLOTO - EFN SOURCE STATEMENT - IFN(S) 



C BOB COFER 3/22/67 

C PLOTS PA *P3 . ANO U COMPONENTS VS. UNIVAR PASSES FCR PATTERN SEARCH 
SUBROUTINE PLOTOP 

DIMENSION ULABEL<10). UPS(500)•XAXIS1500) 

COMMON N.UA(IO).U8{|0).PA.P3*UCIC)*P•PF*UF(10).LC11 
C .B VA(10.K1.K2.K3 »K4•COUNT *STE(10)•MINC.LCl•LC3 

COMMON /PLOT/ US{20»500). NPT .PAS(500) .PBS(500) 

PEAL IN 

DATA (ULABEL(I) •1=1 .10) /6HU(1).. .6HU(2)••.6HU(3).••6HU(4)•••6HUC 
C5}.. .6HU(6). « *6HU( 7) . . .6HU(8).* - 9 6HU { 9} • . .6HU < 10)./ 

C SCALE X AXIS FOR ALL PLOTS 
IN * 0. 

DO 10 I * l.NPT 
IN * IN F I. 

10 XAXISCI) « IN 

5 FORMAT (SE20.8) 

WRITE (6.5) XAXIS.PAS.PBS 

CALL PSCLE (XAXIS.NPT.I0..XMIH.X5CL) 

C PLOT PA AND PB VSc COMPLETE PASSES THROUGH UNIVAR 
C SET INITIAL PA AND PB * Q« TO AVOID EXTREMELY LARGE PB VALUES 
DO U 1 x 1.3 
PASCI) * 0. 

12 PBS(I) * 0. 

CALL PSCLE (PAS.NPT.8..PAMlN.PASCL) 

CALL PSCLE (PBS.NPT.8..PBMIN.P33CL) 

CALL PA^IS (0..0..23H COMPLETED EXPLORATORY MOVES•28 * 10..0.,KMIN. 
CXSCL) 

CALL PAXIS € 0*•0.•1 OH PA AND PB.10•10..90..PSMIN.PBSCL) 

CALL PLINE (XAXIS.PAS.NPT.l) 

CALL PLINE (XAXIS.PBS.MPT.l) 

CALL PLEND (14..0) 

C PLOT U COORDINATES VS. COMPLETE PASSES THROUGH UNIVAR 
DO 15 I * l.N 
DO II J « l.NPT 

11 UPS(J) m US(t.J) 

CALL PSCLE (UPS.NPT.8..UMIN.USCL) 

CALL PAXIS (0..0..28H COMPLETED EXPLORATORY MOVES.28.10.#0•*XWIN# 
CXSCL) 

CALL PAXIS<0..0..ULABEL(I).6.8..90..UMIN.USCL) 

CALL PLINE (XAXIS.UPS.NPT.l) 

15 CALL PLEND <14..C) 

CALL PLEND (l..i) 

RETURN 

END 
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Appendix Y (Unclassified) 

NONLINEAR PROGRAMMING: A REVIEW OF THE STATE-OF-THE-ART 

1. INTRODUCTION 

In recent years, there has been an increasingly widespread recog¬ 
nition that many important practical problems can be formulated as prob¬ 
lems of constrained minimization or maximization. That is, many problems 
of design, control, resource allocation, etc. , can be formulated as follows: 
find the vector argument x which minimizes (or maximizes) a given scalar 
function f(x), subject to a given set of equality and/or inequality constraints. 
The function f is called the objective function. It is also known as the cost 
function, the performance index, the payoff function, etc. 

If the function f is linear, and the constraints are linear equalities 
and/or inequalities, this is the general linear programming problem, for 
which theory and practice are highly developed. If the oojective function or 
the constraints, or both, are nonlinear, it is the general nonlinear program- 
ming problem, also known under various other names: function minimization, 
constrained minimization, mathematical programming, optimization, param¬ 
eter optimization, design optimization, etc. These names have differing con¬ 
notations and associations. n Nonlinear programming' 1 is the term favored 
by mathematicians. "Mathematical programming" is a general expression 
which includes linear, nonlinear, and integer programming, and is princi¬ 
pally used in the field of operations research. The other terms are mostly 
used by engineers. The terms "function minimization" and "parameter opti¬ 
mization" definitely imply that there are only a finite number of variables, 
so the minimization/optimization problem is a problem in the ordinary cal¬ 
culus, as distinguished from the calculus of variations. The other names 
are ambiguous in this respect; they may be used in a broad sense which in¬ 
cludes both finite-dimensional and infinite-dimensional cases, but are more 
often used in a restricted sense which excludes the latter alternative. 
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The development of numerical methods for solution of nonlinear 
programming problems is still in a relatively primitive state, but is a very 
active area of research. Unfortunately, most of the research work has been 
strongly biased toward one or the other of two distinct approaches which have 
not, as yet, been satisfactorily combined,, One of these approaches views 
nonlinear programming as a generalization of linear programming. In linear 
programming the constraints are everything; without constraints there would 
be no problem. Therefore, this approach leads to emphasis on handling 
medium to large numbers of constraints, by modifications of the highly suc¬ 
cessful algorithms for linear programming. The basic assumption made is 
that the objective and constraint functions are approximately linear over 
regions large enough to permit reasonable-size steps toward the solution. 

The other approach starts with consideration of difficult problems in uncon¬ 
strained minimization (e. g, , problems where the objective function has very 
narrow, curving valleys) and adds provisions for constraints as an after¬ 
thought, if at all. An important aspect of the research v/ork based on this 
approach has been the development of methods with seccnd-order (or at least, 
higher than first-order) convergence in the neighborhood of the solution. This 
is not only desirable in itself; it is almost a necessity in cases where the 
matrix of second partial derivatives of the objective function has a wide range 
of eigenvalues. 

The constraint-oriented approach has been the basis of some highly 
successful, large scale optimization programs, such as POP (Ref. 1), COP 
(Ref. 2) and NDS (Ref. 3). These programs have been made highly versatile 
and efficient by sophisticated programming techniques, which may make them 
competitive even in applications where the other approach is more appropriate 
in principle. They are potentially inefficient in cases where very narrow 
valleys occur, but in practical problems of design optimization, process con¬ 
trol, etc. , such cases do not occur very frequently if reasonable care is used 
in problem formulation. Usually, the objective function and constraints are 
quite simple and well-behaved, so constraint-oriented optimization algorithms 
can be highly efficient. However, important exceptions are known to exist. 


AV-3 


By combining the principles of both approaches, it should be pos¬ 
sible to develop algorithms which handle general sets of constraints ef¬ 
ficiently and have accelerated convergence even for unfavorable (but smooth) 
forms of the objective function, or the constraints, or both. Such a com¬ 
bination has not yet been demonstrated in a completely satisfactory form, 
but numerous investigators are working toward this goal, and significant 
progress has already been made. Major advances may therefore be expected 
during the next few years. 

The review briefly describes the nonlinear programming problem 
and some of the principal techniques that have been proposed for its solution. 
Most of the discussion is restricted to the finite-dimensional case, but opti¬ 
mization in function space (i. e. , the optimal control problem) is considered 
in one section. A brief Bibliography is included. It is limited to books, re¬ 
view articles, and technical papers which describe significant technical ad¬ 
vances. 

2. FORMULATIONS OF THE PROBLEM 

The general nonlinear programming problem in n dimensions is to 
find an n-component argument vector x which is feasible (i. e. , satisfies a 
given set of constraints) and minimizes the objective function f(x) subject to 
these constraints. Two extreme special cases are those with no constraints 
(unconstrained minimization) and those where the constraints completely de¬ 
termine x, so the nonlinear programming problem reduces to solving a set 
of nonlinear equations. There are many ways of expressing the constraints. 
One way is given by the relations 


g^x) = 0 

0 < i < m 1 

(la) 

g^x) £ 0 

m^< i < m 

(lb) 


which give m^ equality constraints and (m-m^) inequality constraints. This 
form will be used in most of the discussions below. Other forms may be 
better for computational use. Inequality constraints often occur in pairs, 
giving a limited range for some variable, and this;tail be exploited to reduce 
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the number of separate relations. For example, in the optimization algorithm 
POP (Ref. 1) the constraints are stated in the form 

a^ < x < b (2 a) 

£ < X. < d (2b) 

where a, b, £, d are constant vectors, and y = y (x) is a vector function of 
x, not necessarily of the same dimensionality. Another form, which appears 
different but is equivalent to either of the above formulations, is 

£ < x < b_ . (3a) 

g. (x) = 0 0 < i < m^ (3b) 

This formulation is used in Abadie's Reduced Gradient method (Ref. 4). 

3. NECESSARY CONDITIONS FOR OPTIMALITY 

If f(x) possesses continuous first partial derivatives, a necessary 

condition for an unconstrained local minimum at x is 

—o 

grad f(x ) = 0 (4) 

— o 

The corresponding necessary conditions for constrained cases are 
the Kuhn-Tucker conditions (Ref„ 5) which involve the function 

m 

F (x, X ) = f(x) + 2 X.g.(x) (5) 

1 11 

The Kuhn-Tucker conditions are that the constraint relations given 
by Eq (1) hold, that the gradient of F(x, \ ) with respect to x vanishes, and 
that for values of i corresponding to inequality constraints (m^< i < m) the 
multipliers X^ satisfy 

X. £ 0 if g^(x^) = 0 (6) 

X. = 0 if g. (x ) > 0 

l i — o 

There are also second-order necessary conditions, if the relevant partial 
derivatives exist. For unconstrained minimization, the second-order neces¬ 
sary condition is that M(x), the matrix of second partial derivatives of f(x), 
be positive semidefinite at x^: 

M(x ) > 0 
— o 


( 7 ) 
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The second-order necessary condition for constrained cases is less familiar. 
It is given by Hadley (Ref. 6, page 101-102) in the following form: Let 
(x, \ ) denote the matrix of second partials of F(x, X) with respect to com¬ 
ponents of x, and define G(x) to be the m x n matrix whose ith now is the 
gradient of g^(x) if the ith constraint is active (holds as an equality) and is a 
row of zeros if the ith constraint is inactive (holds as a strict inequality). 

The latter alternative is possible only for inequality constraints.- The neces¬ 
sary condition is that M (x q , X ) d be positive semidefinite for every 

vector d satisfying G(x^)d = 0. An equivalent condition is that some m x m 
matrix N must exist such that 

M* (x X ) = M (x , X ) + G T (x ) NG(x ) (8) 

F —o — F “O — —o —o 

is positive semidefinite. Without loss of generality, N may be assumed, to 

be some multiple of the unit matrix. This necessary condition is related to 

the obvious fact that if x gives a constrained local minimum of f(x), it must 

do the same for 

t” (x) = f(x) + 2 g (x) N g. (x) (9) 

a.c. 1 11 3 

where N is positive semidefinite but otherwise arbitrary, and the summation 

is over constraints active at x . 

—o 

A sufficient condition for a strong local minimum is that the first- 

* 

order conditions hold, and that M (in the unconstrained case) or M (in the 

F 

constrained case) be positive definite, rather than semidefinite only. 

4. THE NWTON- .R APHSON METHOD 

The Newton-Raphson method of solving a nonlinear programming 
problem is based on solving the first-order necessary conditions as a set of 
simultaneous equations and/or inequalities. For unconstrained cases, a 
Taylor expansion of h(x) = grad f (x) gives 

h (x + A x) « h (x) + M(x) A x (10) 

which suggests the Newton-Raphson iteration 

new x = x - M * (x) h (x) 


( 11 ) 
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This iteration converges rapidly in some neighborhood of the minimum. Un¬ 
fortunately, it converges just as readily to stationary points which are not 
minima, and often fails to converge at all. Powell (Ref. 7) has noted that 
convergence can be improved by regarding the Newton-Raphson step as in¬ 
dicating only a direction of change of x, and letting the magnitude of the step 
be determined by a one-dimensional search for the minimum of f(x) along a 
line. However, he found cases in which this also fails, even if f(x) has no 
stationary points other than the minimum. 

If M(x) is positive definite (as it must be, .in the neighborhood of a 
strong local minimum) a sufficiently short move in the direction given by the 
Newton-Raphson correction necessarily decreases f(x). The convergence 
difficulty arises if M(x) has negative eigenvalues for some values of x, en¬ 
abling f(x) to be increased or left unchanged. An obvious remedy is to modi¬ 
fy the inverse matrix in Eq (11), forcing it to become positive definite if not 
already so. An elegant procedure for doing this was invented by Davidon 
(Ref. 8). Davidon's method will be discussed in Section 8 of this Appendix. 

For constrained cases, the procedure is to assume that certain con¬ 
straints are active and others are inactive, ignore the inactive constraints, 
and treat all active constraints as equality constraints. The assumed set of 
active constraints is justified (or corrected) aJ 
pansion of the first-order conditions gives 


0 = 


g 


"M* (x, X ) G T (x)" 


G (x) 


the 

fact. 


1 

> 


AX 


( 12 ) 


where g is a column vector of the constraint functions g.(x), 

2- i — 


h* (x, X ) = h(x) + G (x) 
is the gradient, with respect to x, of 

F* (x, X_ ) = f*(x) + X T g (x) 


X + 


2Ng(x)] 


(13) 


( 14 ) 


and 
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The matrix N is to be chosen (if possible) to make M positive definite. 

F 

Eq (12) suggests a Newton-Raphson iteration analogous to Eq (11). If this 
iteration converges, it converges to a constrained stationary point, which is 

o> 

'I' 

not necessarily a constrained minimum unless can be made positive semi- 
definite by some choice of N. Convergence to non-minimal stationary points 
can presumably be prevented, as in the unconstrained case, by modification 

v»> 

of the matrix M to force positive-definiteness. However, no analog of 
Davidon ! s method, applicable to constrained cases, appears to have been 
published so far. 

5. CLASSIFICATION OF METHODS 

In spite of its attractive feature of rapid convergence near the solu¬ 
tion, the Newton-Raphs on method is generally unsatisfactory because of (a) 
need for second derivatives which may be difficult and/or expensive to com¬ 
pute, (b) convergence uncertainties, including the possibility of false con¬ 
vergence , and (c) in constrained cases, the need to guess which constraints 
are active, with tedious recomputations if the initial guess is wrong c There¬ 
fore, a large variety of other methods have been developed. These may be 
classified according to the way in which they acquire derivative information 
(analytical derivatives, numerical derivatives obtained by finite differences, 
or no explicit computation of derivatives) and according to their speed of con¬ 
vergence near the solution (first order, second order), As with most classi¬ 
fication systems, these classifications are not entirely clear-cut. Some 
algorithms, for example, use analytic first derivatives but compute second 
derivatives by finite differences. Also, many of the algorithms loosely 
described as second order have not been proven to have true second-order 
convergence, although this convergence is faster than first order. The two 
classifications are related, since any method with superlinear convergence 
must make use of second derivatives or equivalent information. For appli¬ 
cations where a large variety of different problems must be solved, algorithms 
which use no derivatives, or compute derivatives by finite differences, are 
generally to be preferred over algorithms which require analytic derivatives, 
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-since the latter require much more input information for specification of a 
new problem. Explicit computation of derivatives by finite differences, as 
compared to implicit use of derivative information from changes of function 
value, has the advantage that the differencing technique can be changed in 
the course of the computation, to provide increased accuracy when needed; 
the algorithm can use simple differences at first, and then go over to 
symmetric differences. 

For applications where essentially the same problem is to be 
solved many times, with slightly different parameter values each time (as in 
on-line optimal control) superlinear convergence is highly desirable not only 
for economy of computation, but because it gives a more smoothly varying 
output. 

6, HEURISTIC DIRECT-SEARCH METHODS 

The two best-known methods which make no use of derivative 
information, and are purely heuristic in motivation, are Rosenbrock's 
method (Reference 9) and the n Direct Search’ or n Pattern Move Search 1 
method of Hooke and Jeeves (Reference 10). Numerous variants of both 
of these methods exist. Neither method, in its original form, included pro¬ 
visions for constraints, but ad hoc modifications for constrained problems 
have been proposed by several authors. 

In Rosenbrock’s method, the program remembers a set of 
orthogonal search directions and uses them in sequence. The relative 
success of the different directions, as measured by net displacements of the 
argument vector, is then used to generate a new set of search directions. 

The procedure is intended to make the search directions tend to align them¬ 
selves with principal axes of the quasi-ellipsoids which are iso-value 
contours of the objective function near a minimum. This behavior has never 
been proved mathematically, but is supported by experimental evidence. 
Convergence is first-order, but usually reasonably rapid. 
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The basic idea of Pattern Move Search is that what succeeds 
once, should be tried again. Starting from an initial ,f base point", an 
exploratory subroutine generates another point with smaller objective 
function. This constitutes a new base point. The vector difference between 
the two most recent base points defines a "pattern move" which is applied to 
the newest base point to generate an "advance point" from which the explor¬ 
atory subroutine is used again. If the overall result of the pattern move 
and subsequent exploration is an improvement, it defines a new base point 
and the process repeats. Otherwise, the procedure restarts from the most 
recent base point, with an exploratory sequence. Pattern Move Search is a 
general framework, into which almost any other algorithm cam be fitted as 
exploratory subroutine. With an appropriate choice for this subroutine, it 
has been found to be an efficient "valley-follower", and gives reasonably 
rapid (but first-order) convergence in most cases. It can handle constraints 
if the exploratory subroutine has this capability, and if the main program is 
modified to enforce feasibility initially and after each pattern move. 

7. GRADIENT METHODS 

For unconstrained cases, the basic equation of the gradient, or 
steepest descent, approach to function minimization is 

new x - x - k H h (x) (15) 

where k is a positive scalar multiplier which may be held constant or itera¬ 
tively adjusted, and the positive-definite matrix H is the "inverse metric" 
i, e. , the matrix whose Liverse defines vector magnitudes by the relation 

(Ax| 2 = Ax^ H * Ax (16) 

The use of a factor like H in Equation (15) is unavoidable since the vectors 
x and h do not transform alike under coordinate transformations, and generally 
do not even have the same dimensionality. The convergence properties of the 
gradient method, can be studied by considering the special case of a quadratic 
function with its minimum at x ~ 0. This function is of the form 



AV -10 


f (x) = 0. 5 x T M q x (17) 

*so its gradient is 

h (x) - x - (18) 

If the vectors w and scalars are the eigenvectors and eigenvalues, res¬ 
pectively, of HM^, x has an expansion of the form 

x = 2 c. v. (19) 

- i i~i 

with scalar coefficients c^. The eigenvectors have the orthogonality property 
v T H v. = 0 ~ v T M v. if i^j (20) 

-l —j -i o -j 

so 

|x| 2 = X |c.| 2 (21) 

Expressing bofh sides of Equation (15) in terms of eigenvectors gives 

new c. = (1 - ka.) c. (22) 

ill 

If all the are nearly equal, i. e. , if HM approximates a multiple of the 
unit matrix, it is possible to choose k so that 

| 1 - kxj « 1 (23) 

for all i, and the gradient method gives rapid convergence. In the more 
usual case where the GL differ greatly from each other (in particular, if the 
magnitudes of the largest and smallest eigenvalues have a large ratio) this 
is not possible, every choice of k either makes some c^ increase in magnitude, 
or leaves other c^ almost unchanged, so a k chosen small enough for stability 
gives very sluggish convergence. Dynamic variation of k from iteration to 
iteration alleviates this situation only slightly. Ultimate convergence of 
f(x) to a local minimum can be assured (if H is positive definite) by con¬ 
trolling k by rules that force monotone decrease of f(x), but computational 
efficiency depends on choosing the metric H to be similar to M q , the matrix 
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of second partial derivatives at the solution-point. Davidon's method, des¬ 
cribed in Section 8 of this Appendix, is an automatic procedure for itera- 

-1 

tively adjusting H to make it approach M (x) if possible, but remain posi¬ 
tive definite always. 

For constrained cases, the gradient method is replaced by the 
"projected gradient 11 method of Rosen (Refs 11, 12) or the "reduced gradient" 
method of Wolfe (Ref 4). In the projected gradient method, the specified 
Ax at each iteration is the sum of two displacements: a displacement which 
attempts to null all deviations from the active constraints, plus another 
displacement parallel to the projection of the gradient vector unto the con¬ 
straint surfaces. In the reduced gradient method, the active constraints 
are used to eliminate some components of x by determining them as 
functions of the remaining components, and the unconstrained gradient 
method is then applied in a smaller-dimensional space. In both these 
techniques, efficient means are provided for treating changes in the set 
of active constraints. However for both methods, efficiency depends on an 
appropriate choice of metric. No automatic procedure for iteratively ad¬ 
justing the metric (analogous to Davidon 1 s method) has yet been proposed 
for constrained cases, although Kelley 1 s method (Ref 13) accomplishes 
essentially the same result. 

Three gradient methods which warrant special discussion have 
already been mentioned in the Introduction. One of these is the NDS pro¬ 
gram of Beskind (Ref 3). This program is based on the projected-gradient 
method. It employs a logarithmic scaling which gives a reasonable metric 
in most practical cases. Also, computation is economized by a method 
of obtaining and updating approximate gradient information for the objective 
and constraint functions, as a byproduct of function evaluations performed 
in the natural course of the search for the optimum. The other two methods 
are POP and COP (Refs 1 and 2) which are slightly different versions of the 
same algorithm, designed for different computers. The algorithm is a 
sophisticated form of the Simplex (LP) method, applied by local lineariza- 
tion. A number of ingenious techniques are employed to minimize computa¬ 
tion. 
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8. SECOND-ORDER METHODS FOR UNCONSTRAINED 

MINIMIZATION 

There are three principal methods of unconstrained minimiza¬ 
tion which are commonly referred to as second-order methods: Davidon's 
method (Ref 8), the Fletcher-Reeves conjugate gradient method (Ref 14), 
and Powell’s conjugate directions method (Ref 15). Strict mathematical 
proofs of second-order convergence for these methods have not been given, 
but their convergence is known to be super-linear (faster than first order). 

In Davidon 1 s method, which has already been referred to, each 
change of x and h (x) is followed by an adjustment of the H matrix, of the 
form 

new H = dj (HAh) (HAh) T + d 2 (Ax) ( Ax) T (24) 

where dj and d^ are computed coefficients. Since Davidon’s paper is not 
easily accessible, his method is most widely known in a slightly modified 
form given by Fletcher and Powell (Ref 16). However, some users of the 
modified method have reported numerical difficulties (loss of positive 
definiteness, due to computational errors) which apparently do not occur 
with Davidon* s original method. 

Both the original and the Fletcher-Powell versions of Davidon’s 
method require analytical evaluation of the gradient vector h(x). A form 
of Davidon’s method using finite differences to determine approximately 
h(x) from changes of f(x) has recently been reported by Stewart (Ref 17). 

This constitutes a significant advance. 

The Fletcher-Reeves method of conjugate gradients uses analytic 
evaluation of the gradient, but could presumably be modified to use finite 
differences instead. Its basic cycle involves 3 vectors: x, h, and s_. 
Initially, s = -h. The cycle steps are 


1 . 

Choose k to minimize f(x.+ ks) 

2. 

New x = x + ks 


3, 

Evaluate new h 

= h (new x) 

4. 

b = |new bp/ | 

hp _ ■> 
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5. New s_ = bs_ - (new h) 

6. Go to 1 

For a quadratic function, this procedure gives the exact minimum in n cycles. 
For non-quadratic functions, it gives rapid convergence, but requires peri¬ 
odic restart every n cycles. A modified version not requiring restart is 
probably possible. 

In Powell 1 s method of conjugate directions, the algorithm remem¬ 
bers a set of n search directions, and performs a one-dimensional minimizing 
search in each of these directions in sequence. The net resulting change of 
x defines a new direction, which (usually) replaces one of the stored directions. 
No analytical derivatives are used, but equivalent information is obtained 
from changes of function value. The method finds the exact minimum of a 
quadratic function in n cycles, and lias rapid (supe'rlinear) convergence for 
non-quadratic functions. However, there are reports of decreased efficiency 
when n is large (above 15 or so). None of these second-mder methods have, 
as yet, been generalized to deal with constraints, except by the use of 
penalty functions. 


9. PENALTY-FUNCTION METHODS FOR CONSTRAINTS 


A convenient and flexible, but approximate, method of enforcing 
constraints in function-minimization problems is by the use of penalty func¬ 
tions. In one version of the penalty-function method, the objective function 
£(x) is replaced by the modified objective function 

f p (x) = f(x) + 2 k. d. 2 (g.(x)) (25) 

where the coefficients K->o are chosen constants, g^(x) is the constraint 
function, and 


<My) 


{ y 'J if the ith constraint f = o' 

min(o,y) j is of the form | g^ > o 


(26) 


Eq (25) could be generalized by using a quadratic form in the 
quantities instead of a weighted sum of squares, but no particular ad¬ 
vantage would be gained by doing so. The quantities d^ 2 in Eq (25) are 
called penalty functions. The larger JL is chosen, the more accurately the 
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i constraint is enforced. Therefore, an approximate solution found with 
a moderate value of can be improved by trying again with a larger value 
of K^. The®displacement from the true solution will be approximately pro¬ 
portional to 1/K^. Therefore, constraints which must be satisfied with 
high accuracy require large values of K.. This is a potential source of 

difficulty, since large values of K. make the function f (x) have "valleys" 

1 P 

with at least one steep wall. The valleys are curved if the constraint- 
surfaces are curved, and not all algorithms for unconstrained minimiza¬ 
tion can follow such valleys efficiently. Algorithms with quadratic con¬ 
vergence (such as Powell's or Davidon's) become highly desirable, though 
other methods (e. g. , Rosenbrock* s or Pattern Move) may also be effective. 
An automatic method for adjusting the coefficients has been given by 
Kelley (Ref 13} who uses penalty functions in conjunction with Davidon's 
method. As a byproduct, Kelley's ingenious method generates an ex¬ 
cellent metric for use with the projectcd-gradient method in a "refinement 
phase" which follows his rough preliminary minimization. 

In the SUMT algorithm (Ref 18) equality constraints are treated 
as described above, but a slightly different approach is used for inequality 
constraints. The penalty function is nonzero even when the constraint is 
satisfied, and varies inversely with distance from the constraint surface. 

On the violation side, the penalty is infinite: moves across this surface 
are rejected as absolute failures. To improve accuracy, the coefficient of 
the penalty function is decreased, rather than increased as for quadratic 
penalties. The advantage of the SUMT method is that each approximate 
solution is feasible, so far as inequality constraints are concerned. 

It is instructive to compare the penalty-function method with 
the Lagrange multiplier method of function minimization. In the latter 
method, we form 

F(x,X) = f(x) + \. gi (x) (27) 

and seek values of the multipliers \., such that the solution of the resulting 
unconstrained problem satisfies the constraints. We also require 
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X^ = 0 for all inequality constraints which are inactive (strongly 
satisfied) 

Xj, - 0 for all inequality constraints of the form g^(x) - 0 
A necessary condition for x Q to give the unconstrained minimum of F(x,X) is 
VF(x 0 ,X) = h(x 0 ) + 2 X i Vgi^) (28) 

Compare this with the result of differentiating Equation (25), which is 
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K i d i(Si 



Vgi(Xo) 


(29) 


so the two are alike if we identify 

2K. d .(g.(x )) = X (30) 

1 I 1 —o 1 


This suggests that if K. is increased, g.~ (x^ and <T (g^ ( Xq )) will decrease 
in magnitude, and the product will approach a limit which is one-half the 
Lagrange multiplier. Therefore, the penalty method may be regarded as 
a servo for adjusting X^ to make it enforce the i tn constraint, and 2IC as 
the feedback gain for this servo. The servo has a standoff error inversely 
proportional to the gain. 

10. COMPARATIVE STUDIES OF NLP ALGORITHMS 

Relatively little has been published about extensive, objective 
comparisons of different algorithms for nonlinear programming. Numerous 
authors have reported comparisons between their proposed new methods 
and one or more previously existing methods, with invariably favorable con¬ 
clusions. A number of reasonably objective comparisons of selected 
methods have been reported (see review articles in the Bibliography). An 
extensive survey has been conducted by Colville (Ref 19) by sending a 
questionnaire and set of problems to possessors of optimization algorithms, 
with requests for 

(a) additional test problems 

(b) description of algorithm 

(c) performance data for the test problems (success or not, 
preparation time, .running time, > number of evaluations, etc.) 
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Several dozen replies were received. To correct for dif¬ 
ferences in computer speed, each reported running time was normalized 
by dividing by the time required for execution, on the same computer, of a 
standard computation defined by Colville. However, even the normalized 
times are difficult to interpret, since no standard was given for the degree 
of accuracy to be sought. More meaningful, perhaps, is the fact that for 
most of the methods, success was not reported for all problems. However, 
this may reflect limited effort rather than limited capability. Since all but 
one of Colville*s problems involve constraints, algorithms designed purely 
for unconstrained minimization rated poorly in his comparisons (and rightly 
so). 

11, NONLINEAR PROGRAMMING FOR TRAJECTORY 

OPTIMIZAT ION 

There are two distinct ways of applying nonlinear programming 
to problems of optimal control of dynamic systems (trajectory optimization). 
The first way, which is called the n indirect M method for historical reasons, 
is to make use of necessary conditions for optimality derived from the cal¬ 
culus of variations (or its modern counterpart, optimal control theory) to 
reduce the candidate trajectories to an n-parameter family of extremals 
defined by given sets of differential equations. The objective function and 
the constraint functions then become ordinary functions of n parameters, 
and the algorithms for finite-dimensional NLP are applicable. Each evalua¬ 
tion of the objective function or the constraints requires integration of a 
set of differential equations. 

The second way, which is called the "direct" method, is to use 
NLP in function space. The n-component argument vector x is replaced by 
a control history u (t), regarded as a vector u in function-space, and the 
objective function becomes a functional of u. The constraints may involve 
instantaneous values of u(t), or functionals of u, or both. Many of the NLP 
algorithms for n dimensions have analogs in function space. For example, 
the function-space analog of the n-dimensional Newton-Raphson method'has 
been popularized byBellmannas his n quasilinearization n method, and by 
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McGill and Roberts as the n Generalized Newton-Raphson" method (Ref 20). 

The fT Steepest ascent" methods of Bryson and Kelley (Refs 21, 22) are 
analogous to the projected gradient method in n variables. As might be ex¬ 
pected, these methods are critically dependent on the choice of metric. 

The convergence properties can be analyzed by considering the case where 
the objective function is a quadratic functional of u, written symbolically as 
£ = 0. 5 u T Mu (31) 

where M is a linear operator in function space. The inverse metric is 
defined by some other linea.r operator H, which must have an easily com¬ 
puted inverse. Then 

|u| 2 = u T H~ k (32) 

The convergence analysis sketched in Section 7 can be carried out in func¬ 
tion space by considering the eigenvalue problem 
CtH * v = Mv 

There will generally be an infinite number of eigenvalues. If their magni¬ 
tudes vary over a wide (or infinite) range, convergence will be poor (or im¬ 
possible). Under fairly general conditions, an especially favorable choice of 
the metric is 

r p I A t <~p 

u H u = J ^ u (t)W(t)u(t)dt (33) 

u 

where W(t) is the matrix of partial derivatives of the variational Hamiltonian 
with respect to control variables. With this choice of metric, it can be 
shown that the eigenvalues approach a finite limit* so all but a finite number 
of them will be within a small range. This greatly improves the convergence 
properties, and is the theoretical justification of the n Min-H n method of 
Kelley (Ref 22) and Gottleib (Ref 23). In favorable cases, the n Min-H n 
method can give excellent convergence. Gottleib has generalized the "Min-H 1 ' 
method to cases which cannot be handled by the metric given by Eq, (33). His 
generalization involves certain matrices which must be chosen arbitrarily, and 
critically affect the convergence. No general rules are available for selecting 
these matrices, but Gottleib has been able to find good choices for his example 


cases. 
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A function-space analog of the Fletcher-Reeves conjugate gra¬ 
dient method has been proposed by Mitter and Lasdon (Ref. 24). The second 
order convergence of the conjugate gradient method in n dimensions does not 
carry over into function space, but a major improvement over ordinary gra¬ 
dient methods can be expected. Theory indicates that this method, modified 
by use of the n Min~H n metric, should be capable of convergence rates not 
greatly inferior to a second order method. 
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